Advection-dominated transport past isolated disordered sinks: stepping beyond homogenization

We investigate the transport of a solute past isolated sinks in a bounded domain when advection is dominant over diffusion, evaluating the effectiveness of homogenization approximations when sinks are distributed uniformly randomly in space. Corrections to such approximations can be non-local, non-smooth and non-Gaussian, depending on the physical parameters (a Péclet number Pe, assumed large, and a Damköhler number Da) and the compactness of the sinks. In one spatial dimension, solute distributions develop a staircase structure for large Pe, with corrections being better described with credible intervals than with traditional moments. In two and three dimensions, solute distributions are near-singular at each sink (and regularized by sink size), but their moments can be smooth as a result of ensemble averaging over variable sink locations. We approximate corrections to a homogenization approximation using a moment-expansion method, replacing the Green’s function by its free-space form, and test predictions against simulation. We show how, in two or three dimensions, the leading-order impact of disorder can be captured in a homogenization approximation for the ensemble mean concentration through a modification to Da that grows with diminishing sink size.

We investigate the transport of a solute past isolated sinks in a bounded domain when advection is dominant over diffusion, evaluating the effectiveness of homogenization approximations when sinks are distributed uniformly randomly in space. Corrections to such approximations can be non-local, nonsmooth and non-Gaussian, depending on the physical parameters (a Péclet number Pe, assumed large, and a Damköhler number Da) and the compactness of the sinks. In one spatial dimension, solute distributions develop a staircase structure for large Pe, with corrections being better described with credible intervals than with traditional moments. In two and three dimensions, solute distributions are near-singular at each sink (and regularized by sink size), but their moments can be smooth as a result of ensemble averaging over variable sink locations. We approximate corrections to a homogenization approximation using a moment-expansion method, replacing the Green's function by its free-space form, and test predictions against simulation. We show how, in two or three dimensions, the leading-order impact of disorder can be captured in a homogenization approximation for the ensemble mean concentration through a modification to Da that grows with diminishing sink size.

Introduction
Transport processes in many natural systems take place in spatially disordered domains. In many instances, these processes can be adequately described by averaging procedures, Darcy's law describing flow in random porous media being a well-known example [1]. However, it is important to understand the impact of disorder, particularly in instances where disorder has a significant influence (for example, in explaining breakthrough effects, whereby solute is carried rapidly along a small number of high-flow paths through a random porous medium [2]). The present study contributes to this effort by characterizing the impact of spatial disorder on the uptake of a solute that is advected past distributions of isolated sinks. This problem is loosely motivated by transport of maternal blood in the intervillous space of the human placenta [3] but is posed here in more general terms.
A common assumption that is exploited in order to describe transport in media with complex microstructure is to assume periodicity at the microscale [4][5][6][7]. This allows an asymptotic twoscale expansion to be developed, with a unit-cell problem (with periodic boundary conditions) being solved in order to provide a description of slowly varying (homogenized) variables at the macroscale. While this approach has been extended to accommodate slow spatial variation of the microscale field [8][9][10] and developed for a variety of applications [11][12][13][14][15][16], it is less adaptable to situations where the microscale exhibits appreciable spatial disorder. Approaches currently adopted in such instances include formal methods of stochastic homogenization [17], spatial averaging techniques [18] or simulations using random microstructures realized within periodic unit cells [19].
A spatially disordered medium can be characterized as a random field with prescribed statistical properties. The 'forward' problem that we address here seeks to understand how these properties map to the statistical properties of the concentration field of a solute as it passes through the medium. This map is mediated by physical processes embodied in a partial differential equation (in the present instance, a linear advection-diffusion-reaction equation). The primary question addressed by a homogenization approximation is how to translate the first moment of the sink density to the first moment of the associated concentration field (where first moments are ensemble averages). More refined questions address the impact of spatial disorder, captured in the second moment (covariance) of the sink density, on the mean and covariance of the concentration field. Provided solute fluctuations are bounded in an appropriate sense, these corrections can be evaluated by perturbation around the leading-order homogenization approximation, as we illustrate below, and as demonstrated previously by Dagan [20], Cushman et al. [21], Chernyavsky et al. [3], Russell et al. [22] and Russell & Jensen [23]. If fluctuations become sufficiently large, or if distributions become strongly non-Gaussian, higher moments (or even full probability distributions) of the solute field may need to be evaluated.
Homogenization approximations exploit the separation of lengthscales between the microscale and the macroscale. However, when considering solute uptake at isolated sinks, a further lengthscale needs consideration. The microscale involves two lengthscales, an intersink distance ρ (assumed small compared with the overall size of the domain) and a sink size ς . As ς becomes vanishingly small with respect to ρ, over the shortest lengthscales, diffusion can be expected to dominate advection in the neighbourhood of sinks, and the concentration field can be expected to be described locally by the solution of a diffusion equation in the neighbourhood of a point source. In one dimension, this leads to a concentration field with a staircase structure, with a thin diffusive boundary layer forming upstream of each sink [23]. In two and three spatial dimensions, large solute gradients surround the sink, and the concentration field grows in magnitude proportionally to log(ρ/ς) and ρ/ς, respectively. This effect amplifies fluctuations, as we demonstrate below, and is known to restrict the applicability of homogenization approximations in two and three dimensions [14].
The present study develops an approach initiated by Russell & Jensen [23], who used an iterative method to approximate the effects of disorder in a linear transport problem involving advection, diffusion and solute uptake via first-order kinetics. They considered a spatially one-dimensional problem with uptake taking place at isolated point sinks. They considered parameter ranges for which a steady concentration field can be constructed via a smooth (homogenized) leading-order solution, to which corrections are added that account for the discreteness and disorder of the sink distribution. Corrections are non-local and were evaluated using a Green's function, sidestepping the assumption of unit-cell periodicity that underlies traditional two-scale homogenization. Russell & Jensen [23] considered a parameter regime in which diffusion was dominant at the intersink distance ρ, allowing the use of Riemann sums to approximate certain sums as integrals. Their approach was constructive: rather than seeking to prove formal convergence, explicit evaluation of the magnitude of corrections allowed domains of validity to be established, and simulation was used to evaluate accuracy. Russell & Jensen [23] demonstrated improved accuracy of corrections to a leading-order homogenization solution evaluated using a Green's function approach in comparison with classical two-scale asymptotics assuming microscale periodicity. They also compared the magnitude of corrections with solute fields for periodic, normally perturbed and uniformly random sink distributions, each showing distinct dependence on the underlying physical parameters.
Here, we extend this work in four directions, while adopting the same constructive approach: (i) the problem is reformulated to focus on the mapping from statistical moments of the sink distribution to statistical moments of the solute distribution, allowing sink distributions to be represented (for example) as a Gaussian process; (ii) a parameter regime is considered for which advection dominates diffusion over intersink lengthscales, leading to non-smooth concentration profiles; (iii) the study is extended to two and three dimensions, for which the point-sink approximation must be relaxed to allow sinks to have finite size, so that fluctuations remain bounded; (iv) although corrections to a naive homogenization approximation are generally nonlocal, we show that an essentially local correction to the mean concentration field can be identified when the sink correlation length is sufficiently small, and we evaluate this correction explicitly for sinks distributed uniformly randomly in a two-or three-dimensional domain.
To set the scene, figure 1 shows a set of realizations of a one-dimensional advection-uptake process (with no solute diffusion). In this example, 19 point sinks are distributed randomly in the domain (0, 1), each removing a fixed proportion of the oncoming concentration (which takes the value 1 at the inlet at x = 0 and is swept uniformly in the positive x direction). An individual realization (magenta) reveals the staircase structure of a typical one-dimensional concentration field and shows how it deviates appreciably from the discontinuous sample median (green) and the smooth sample mean (red). This example illustrates how the concentration distribution can be non-Gaussian, with credible intervals (cyan) deviating from the equivalent intervals defined by the sample variance (blue) near the source (where concentrations cannot exceed unity) and near the sink (where concentrations cannot fall below 1.05 −19 ≈ 0.396). This example illustrates how averaging leads to non-smooth concentration fields having smooth statistical moments, even if these must be interpreted cautiously in some circumstances. Expressions for the moments and credible intervals of this simple example are derived in appendix A.
While it is relatively straightforward to make use of an exact Green's function for a onedimensional transport problem (satisfying appropriate inlet and outlet boundary conditions), this is less true in two and three dimensions, and the high-dimensional integrals needed to evaluate higher moments quickly become computationally costly. However, when advection dominates diffusion, the free-space Green's function provides a potentially useful simplification. The Green's function for advection/diffusion/uptake has a discontinuity in one dimension, a log r singularity in two dimensions and a 1/r singularity in three dimensions, making homogenization feasible for point sinks in one dimension [24] but more challenging in higher dimensions [14]. Accordingly, we consider below isolated sinks of finite width ς , taking them to be distributed uniformly randomly in space. We formulate a transport problem in a domain that is bounded in the advective direction x 1 , assuming a spatially uniform inlet flux at x 1 = 0, and assume that sink distributions are statistically uniform over a region that is bounded in the transverse direction. Despite individual realizations having a complex spatial structure, moments typically depend on x 1 alone, and become smooth as a result of averaging. In the present study, we assume that advection is uniform, ignoring heterogeneity of the flow field or of diffusivity, allowing us to exploit a tractable free-space Green's function.
In order to capture the effect of disorder within a homogenization approximation, we also adopt a device described by Noetinger et al. [25] and exploit the limit in which the correlation length of the covariance of the sink distribution is very small. In the present example, we show that this length is provided by the sink size ς for sinks distributed uniformly randomly in two or three dimensions. This allows us to evaluate an effective uptake parameter Da eff : replacing the dimensionless Damköhler number in the naive homogenized solution with Da eff , we obtain a direct approximation for the mean concentration that quantifies how disorder reduces uptake when sinks are distributed uniformly randomly in two or three dimensions.
The model that we investigate is outlined in §2a, with example simulations presented in §2b. The moments-based expansion is presented in §2c, revealing the critical roles of the Green's function ( §2d) and its singularities in the evaluation of high-dimensional integrals ( §2e). The derivation of Da eff is given in §2f. Predictions are evaluated against simulations in §3.

Model and methods (a) The model problem
We formulate the model in three dimensions, adopting analogues in one and two dimensions when required. Let D 3 be a domain of thickness L defined such that , U, D and S represent the (dimensional) solute concentration field, uniform advective velocity in the x * 1 direction, diffusion coefficient and uptake rate, respectively. Uptake is mediated by a distributed sink function satisfying 1 +ĝ * (x * ; ω) ≥ 0, whereĝ * has zero spatial average. ω denotes thatĝ * (x * ; ω) is a realization drawn from a prescribed distribution, making C * (x * ; ω) a random variable.
We prescribe a solute flux q on the plane x * 1 = 0, with zero diffusive flux on x * 1 = L and as x * 2 , x * 3 → ±∞. Defining x = x * /L,ĝ(x; ω) =ĝ * (x * ; ω) and C(x; ω) = C * (x * ; ω)/(q/U), the dimensionless concentration satisfies the advection-diffusion-uptake equation and boundary conditions . The Péclet number Pe = UL/D represents the strength of advection to diffusion; the Damköhler number Da = SL 2 /D relates the rate of uptake to diffusion. We focus here on the strong-advection regime Pe max(1, √ Da); of particular interest is the distinguished limit in which Pe / Da = U/SL = O(1), implying a balance between advection and uptake across the whole domain.
Isolated sinks are taken to be of finite size and to occupy a subdomain D s 3 of D 3 in which x 1 ∈ [0, 1] and x 2 , x 3 ∈ [−L s , L s ]. Let ρ = 1/N be the average inter-sink distance in any direction, where N ∈ Z + represents the number of sinks per unit length. Let the midpoint of sink locations be represented by Thus there are (2M + 1) 2 /ρ sinks in D s 3 with an average density per unit volume given by ρ −3 . We defineĝ(x; ω) to bê has a spatially averaged density of zero within D s 3 . We assume throughout that isolated sinks have multivariate uniform distribution, such that We adopt the Gaussian sink structure function where ς remains sufficiently small to satisfy (2.3) and prevent sinks from overlapping, to exponential accuracy. This function is chosen for convenience but could be replaced to model specific applications. It will be helpful to represent distributions of isolated sinks in terms of their first two statistical moments. As shown in appendix B, uniformly random sinks with Gaussian structure function (2.4) have ensemble mean and covariance and K represents covariance. An important distinction between one-dimensional and higher-dimensional cases is evident. For n = 1, Kĝ has a non-local contribution (with N sinks in a one-dimensional domain, finding one sink at a location reduces slightly the chance of finding another elsewhere). However, for n > 1, with M → ∞, the non-local term vanishes (because the sinks can occupy an arbitrarily wide area or volume within D 2 or D 3 ). The sink density in this case resembles a Gaussian process with square-exponential covariance σ 2 exp(−|x − y| 2 / 2 ), having variance and correlation length given respectively by   arises through a combination of averaging effects and strong advection, which limits the degree of lateral diffusive spread downstream of each sink. We seek approximations of these smooth onedimensional functions in terms of the sink density ρ, sink width ς and the physical parameters Pe and Da.

(c) A moments-based expansion
The volume-or ensemble-averaged sink density in (2.1a) is unity, making it natural to propose the leading-order homogenized linear and boundary operators associated with (2.1) as respectively. The leading-order homogenized solution C H (x) associated with (2.1) can be found by solving It is evident that C H (x) depends only on x 1 , being where φ ≡ Pe 2 /4 + Da and ψ(x 1 ) ≡ (2 Pe φ + Pe 2 +2 Da) e φx 1 + (2 Pe φ − Pe 2 −2 Da) e −φx 1 . In showing how the concentration decays over a lengthscale defined by a balance between uptake and advection. Writing the concentration as we construct a solution of (2.1), to be validated a posteriori, using the ansatz and etc. To invert the linear operators in (2.10), we define G 3 (x, x ) to be the associated threedimensional Green's function satisfying Applying homogeneous boundary conditions in the x 2 -and x 3 -directions is appropriate as the source term is compact. The Green's function can then be used to give the corrections and We characterize the corrections in terms of their moments evaluated over realizations, specifically and This approach extends to n = 1, 2 dimensions, replacing D 3 and G 3 (x, x ) with D n and G n (x, x ), respectively, generalizing the one-dimensional formulation in Russell & Jensen [23]. In higher dimensions, complications emerge due to singularities of G 2 and G 3 as x → x and the high dimensionality of the quadrature.

(d) The free-space Green's function
While the Green's function in one dimension is straightforward to evaluate (appendix C), it is convenient to instead use the free-space Green's function G n (x − x ) for computations in higher dimensions. In three dimensions, this satisfies G n is given by (C 3): it shares with G n the log(φ|x − x |) singularity in two dimensions and 1/|x − x | singularity in three dimensions. G n offers a close approximation of G n in the limit Pe max(1, √ Da), as illustrated for n = 1 in figure 3a,b. This shows a discrepancy between G 1 (x 1 , x 1 ) and G 1 (x 1 − x 1 ) only within a 1/ Pe distance of the outlet in x 1 and the inlet in will allow us to make use of G 1 later on. G(x − x ) denotes the field in the x plane generated by a point sink at x . In two dimensions [three dimensions], concentration contours have an approximately elliptical [ellipsoidal] shape, with dimensions illustrated in figure 3c when Pe max(1, √ Da), as explained in appendix C. We can use this structure to identify the asymptotic region of influence associated with a point x, within which sources at x will contribute appreciably to the concentration field at x, as illustrated in figure 3d. Strong advection implies that the region of influence is largely upstream of x, while strong uptake ensures that the region is narrow in the direction transverse to the flow. This allows quadrature to be restricted to physically relevant domains. x 1 x 2 x 1

(e) Evaluation of moments
Adopting the free-space Green's function approximation and incorporating the sink moments (2.5), (2.13) becomes and (2.15c) We now consider approximations when the domain width is large (L s ρ) and the sink width small (ς → 0). To approximate the variance of C 1 in this limit, we can replace F (2.14), giving (2.16) This reduces the 2n-dimensional integral (2.15b) to a cheaper n-dimensional integral (2.16), although some loss of accuracy is anticipated by ignoring the finite sink size. While (2.14) can also be used to reduce the second integral in (2.15c) to one dimension, a δ-function approximation cannot be used for the first integral in E[ C 2 ] because of singularities in G 2 and G 3 . Instead, we exploit the fact that F can be approximated by its leading-order singular form. We summarize the results of this calculation (see appendix D), as ς → 0 in n dimensions, as and γ is the Euler-Mascheroni constant. The correction in one dimension is independent of the sink size ς as ς → 0, whereas in two and three dimensions the correction grows in magnitude as ς becomes asymptotically small. In two and three dimensions, when L s ρ, the final terms of In two dimensions, downstream influence can again be captured approximately by using the farfield approximation (C 7) of G 2 in the first integrals of (2.16) and (2.17a) (taking Pe max(1, √ Da), ς 1/ Pe and M → ∞) to give In three dimensions, the same approach using (C 10) yields

(f) Defining the effective Damköhler number
In addition to calculating the mean correction directly via (2.15c), we consider how the homogenization problem can be adjusted to capture the leading-order effect of disorder. We seek the constant Da eff such that the solution of approximates E[C(x; ω)] to a suitable degree of accuracy. The exact solution of (2.21) is identical to the leading-order homogenized solution given in (2.8) but with Da replaced with Da eff , namely Assuming the correction C(x; ω) is small compared with C H , the linear operator can be inverted to give where the ω notation is dropped as the leading-order correction is deterministic. We then rewrite (2.15c) as Comparing this with (2.23) gives the approximate relation . In (2.25), we have expanded the domain D 3 to R 3 , a reasonable assumption when sufficiently far from boundaries and the decay lengthscale of G 3 Kĝ is sufficiently short. Exploiting the fact that Kĝ(x, y) = Kĝ(x − y) depends on x − x rather than x and x independently, we can rewrite (2.25) as where * denotes convolution. If the decay lengthscale in Kĝ is sufficiently short, then G 3 Kĝ resembles a δ-function with the appropriate weight and is given by [25]  For a sink covariance function of the form σ 2 exp(−|x − y| 2 / 2 ), taking → 0 and accounting for the singularity in G 2 and G 3 , we evaluate (2.28) using methods given in appendix E to give where we have included the corresponding one-dimensional approximation using (2.28). Recall that φ = Pe 2 /4 + Da. The correction to Da in (2.29) is proportional to (one-dimensional), 2 log (two-dimensional) and 2 (three-dimensional), showing how the difference between Da and Da eff decreases with dimension for fixed variance and fixed correlation length. In one and two dimensions the correction is proportional to 1/φ and log(φ), respectively, whereas in three dimensions φ does not appear in the correction, demonstrating how the impact of advection on the effective uptake decreases as the dimension size increases. For uniformly random sinks in two and three dimensions letting L s → ∞, we can now use (2.6), noting that the variance depends on sink size, to obtain

Results
The variance of the concentration field in one and two dimensions is illustrated in figure 4. The variance is smooth in both cases, due to strong mixing of sink locations over realizations. In one dimension, because exactly N sinks are encountered along the domain, the concentration at the outlet is strongly constrained (as it was in figure 1), and the variance falls close to zero at the outlet. In two dimensions, this constraint is weaker (N sinks are encountered on average between x 1 = 0 and x 1 = 1), so that the variance remains large at the outlet; (2.19a), for example, predicts that the two-dimensional variance is largest at the outlet for 4 Da < Pe. The cloud plot in figure 4b demonstrates the magnitude of the sampling error from 10 4 two-dimensional simulations, and the independence of the transverse coordinate x 2 (figure 2c). Figure 4a shows how the variance in one dimension predicted by (2.13b) matches closely with the sample variance taken from Monte Carlo simulations. In one dimension, the limit ς → 0 can be taken straightforwardly, using (2.16), and it provides a good approximation to the sample variance and the full integral (2.15b), while overpredicting the predicted variance uniformly. The approximation (2.18a), using the leading-order approximation of the free-space Green's function for Pe 1, captures the shape of the variance well but overestimating its maximum (predicting 0.0081 at x 1 ≈ 0.38 for the chosen parameter values, capturing its x 1 -location well but overestimating its value 0.0063 by almost 30%). In two dimensions, numerical evaluation of (2.15b) is expensive so we show only the simplified approximation (2.16), which overestimates the sample variance by approximately 10% (due to neglect of finite sink size) but captures the overall features reasonably well. The cruder prediction (2.19a  compensates for the leading-order homogenized solution over-predicting uptake. The corrections grow with dimension, particularly through the factors β n from (2.17b) as ς → 0. In two and three dimensions when taking the limit M ρ (i.e. L s is asymptotically large), E[ C 2 (x; ω)] can be simplified as the second integral becomes asymptotically small. Therefore, the computational expense of calculating the correction is further reduced to solving one simple one-dimensional integral. The simpler estimate (2.18b) places the maximum one-dimensional correction within the domain (but downstream of the maximum variance), of O(ρ Da 2 / Pe 2 ). The two-and threedimensional estimates (2.19b,c) place the maximum correction within the domain for Da > Pe, but at the outlet otherwise (as in figure 5), although they do not capture the weak boundary layer near The mean concentration in two and three dimensions predicted using the Da eff approximation (2.30) is shown in figure 5b. The correction to C H grows with dimension, as expected, due to the increasingly large concentration fluctuations near each (regularized) sink. The approximation provides close agreement to Monte Carlo sampling in two dimensions, and to the prediction (2.17a) in two and three dimensions. (Monte Carlo simulations in three dimensions were not attempted.)

Discussion
This study has characterized the impact of spatial disorder on a transport process described by a linear advection/diffusion/uptake equation, assuming a uniformly random distribution of isolated sinks. Using a leading-order homogenization approximation (2.7), (2.8) as a starting point, corrections were computed that describe the likely size of solute fluctuations around a mean field in a particular realization, and the correction due to disorder when evaluating the ensemble mean concentration. Bearing in mind the limitations of using statistical moments to characterize non-Gaussian concentration fields (figure 1), we used a moments-based expansion to relate the mean and covariance of the sink distribution to the mean and covariance of the solute field (2.13). The first two moments of the sink distribution, when sinks are distributed uniformly randomly (2.5), show an important distinction between one dimension and higher dimensions, namely that in a sufficiently wide domain in two and three dimensions the correlation length of the sink covariance is set by the sink width (2.6). Simulations (figure 2) reveal the multiscale nature of the problem: despite large concentration fluctuations in the neighbourhood of individual sinks in an individual realization, ensemble averaging leads to smooth moments of the solute distribution with primary dependence only on a single spatial coordinate. Nevertheless, moments demand calculation of high-dimensional integrals, which we simplified by replacing the exact Green's function with its (explicit) free-space form, confining quadrature to appropriate regions of influence (figure 3d), replacing the regularized sink distribution (where possible) with its δ-function approximation, and integrating over lateral dimensions using identities such as (2.14). This allowed accurate predictions of concentration means (figure 5a) in one and two dimensions, and of variance in one dimension (figure 4); the over-prediction of solute variance in two dimensions would likely be corrected by use of the regularized sink distribution, albeit using more expensive quadrature. Cruder analytical estimates (2.18), (2.19), (2.20) were achieved by neglecting any upstream influence of one sink on another.
For vanishingly small sinks (the limit ς → 0), concentration fields are discontinuous in one dimension (figure 1), and have log(1/ς ) and 1/ς singularities in two and three dimensions, respectively. These appear both in corrections to the ensemble-averaged mean concentration (2.17) and in the effective Damköhler number (2.30) that can be used in a modified homogenization approximation in two and three dimensions. The latter approximation cannot be applied for uniformly random sinks in one dimension, because the sink locations are correlated over the whole domain; however it can be applied when sinks are described by a Gaussian process with sufficiently short correlation length (2.29). For sink distributions of fixed variance, the impact of disorder falls as the sink correlation length vanishes (2.29); however, for uniformly random sinks in two and three dimensions the variance of the equivalent Gaussian process rises as ς falls (2.6), contributing to the reduction in uptake captured in (2.30). As reported by Russell & Jensen [23] and Price [27], mean correctors derived assuming periodic sink distributions show different dependence on parameters to those reported in (2.18b)-(2.20b). For example, considering the expansion (2.9), the dominant corrector in the deterministic periodic problem appears at O(Da) and shares the wavelength of the microstructure, whereas the dominant correction in the uniformly random case is stochastic with smooth variance ((2.18a)-(2.20a), figure 4) with the mean correction appearing at O(Da 2 ) ( figure 5).
This study has a number of obvious extensions, prominent among which is consideration of other types of spatial disorder. For flow in porous media, one expects the flow field to have disorder that correlates appreciably with the disorder in the sink distribution [28]. The present perturbative approach provides a route for understanding the contributions of flow, sinks and their combination to solute distributions, and it will be interesting to evaluate the present approach against predictions of existing studies of reactive transport in porous media relying either on periodicity assumptions [4,6,7] or averaging procedures [18]. Other obvious factors to consider include unsteady effects, variability in sink strength (considered in one dimension by [22]) and the impact of a nonlinear uptake kinetics (considered by Dalwadi & King [12] using two-scale homogenization). As demonstrated by Chernyavsky et al. [3] and others, the statistical properties of the underlying spatial disorder interact with the physical lengthscales associated with transport processes to give a range of possible outcomes. The present study illustrates some of the challenges in stepping away from traditional two-scale approaches towards non-local calculations, drawing attention to the need for efficient schemes for high-dimensional quadrature in order to characterize uncertainty.
We can revisit the expansion (2.9) and use evidence that terms Da C 1 or Da 2 C 2 become comparable in magnitude with C H as evidence of the breakdown of a homogenization approximation. In one dimension, based on the estimates in (2.18), the restriction Pe max(1, √ Da) must be extended to Pe max(1, √ Da, ρ 1/2 Da), which holds along the distinguished limit Pe ∼ Da for arbitrarily large Pe. The parameter Da 2 ρ/ Pe, measuring the magnitude (relative to C H ) of the fluctuation variance and the correction to the mean, takes the value 0.05 in figure 1 (with Pe → ∞, but with Da ρ 1/2 / Pe = S 1 /ρ 1/2 ; see appendix A) and figures 4a and 5a. In these examples, fluctuations with standard deviation of order 20% dominate the correction to the mean, of order 5%. In two and three dimensions, however, the range of validity of the approximation is reduced and the correction to the mean (that grows with diminishing sink size) overtakes the fluctuations as the dominant correction. In three dimensions, we require Pe max(1, √ Da, Da 2 ρ 3 /ς ) (for ρ 3 ς ρ 1), which confines the distinguished limit to 1 Pe ∼ Da ς/ρ 3 . The example shown in figure 5 has ς/ρ 3 = 1.25: as the figure indicates, the predicted correction to the mean is sufficiently large to call into question the validity of the homogenization approximation in this case. In two dimensions, the constraint on the distinguished limit is 1 Pe ∼ Da 1/(ρ 2 log(ρ 2 /ς )): the example in figures 2, 4b and 5a with ρ −2 = 25 therefore sits at this upper threshold, although the predicted corrections are still effective.
All authors gave final approval for publication and agreed to be held accountable for the work performed therein.
Competing interests. We declare we have no competing interests. Appendix A. One-dimensional concentration profiles with zero diffusion Let ξ i (i = 1, . . . , N) denote point sink locations, distributed as order statistics U j:N taken from a uniform distribution U ∼ U (0, 1) with probability density function (pdf) π U (x) = 1 for 0 ≤ x ≤ 1 and zero otherwise. Each sink location follows a Beta distribution [29] such that ξ j ∼ β(j, N − j + 1), where j = 1, . . . , N. Here β(x, y) ≡ t x−1 (1 − t) y−1 /B(x, y), where B(x, y) ≡ Γ (x)Γ (y)/Γ (x + y). The cumulative distribution function (cdf) F ξ j (x) = P(ξ j ≤ x) is given by the regularized incomplete beta function The one-dimensional concentration distribution that falls by a factor 1/(1 + S 1 ) at each sink from its inlet value C 0 = 1 (figure 1) satisfies (This problem can be defined as a limit of the one-dimensional form of (2.1), with Pe → ∞ taking S 1 = Da ρ/ Pe with Da / Pe = O(1).) The probability of being at concentration C j for some given and Var[C(x)] are plotted in figure 1 using (A 1). The cdf of the concentration C j is given by N).
Let the cdf take a value F C j (C) = r. Then (A 6) can be inverted to give the corresponding sink locations asξ N).
We can therefore use (A 2) to find the cdf credible intervals as Credible intervals that ensure 95% of concentration profiles are contained between the two bounds are shown in figure 1 using r = 0.025 and r = 0.975 in (A 8); the median is evaluated using r = 0.5. Credible intervals respect the requirement that the concentration is bounded between C N at the outlet and C 0 at the inlet, demonstrating that the solute distribution is non-Gaussian.

( C 1 )
It is convenient to re-express this to allow numerical evaluation when Pe 2 Da. We expand exponential terms to obtain From the first term inG − (G + ), we see a boundary layer of width approximately 1/ Pe (Da / Pe) exists upstream (downstream) of x 1 = x 1 (figure 3a). The final two terms inG − andG + account for the boundary conditions, which gives a boundary layer of width approximately 1/ Pe at the x 1 -outlet and x 1 -inlet. The n-dimensional free-space Green's function G n (x − x ) associated with (2.1) satisfies L n G n = δ(x − x ) and decays in the far field. Seeking a solution of the form G n (x) = e λx 1 f (r) where r = |x| and setting λ = Pe /2 leads to where φ ≡ Pe 2 /4 + Da and K ν represents the modified Bessel function of the second kind [30]. The free-space Green's function in one dimension is readily evaluated, noting that K ±1/2 (z) = π/(2z) exp(−z), as for suitably small . Using (D 1), we obtain G 2 K(0) ≈ − 1 4 σ 2 2 (γ − 2 log(φ )), hence yielding the two-dimensional result in (2.29). The analogous integral in three dimensions reduces to Again for small , we use (E 2) to evaluate the integral for r = O( ), leading to G 3 K(0) ≈ − 1 2 σ 2 2 , the three-dimensional limit in (2.29).