Stochastic transport in the presence of spatial disorder: fluctuation-induced corrections to homogenization

Motivated by uncertainty quantification in natural transport systems, we investigate an individual-based transport process involving particles undergoing a random walk along a line of point sinks whose strengths are themselves independent random variables. We assume particles are removed from the system via first-order kinetics. We analyse the system using a hierarchy of approaches when the sinks are sparsely distributed, including a stochastic homogenization approximation that yields explicit predictions for the extrinsic disorder in the stationary state due to sink strength fluctuations. The extrinsic noise induces long-range spatial correlations in the particle concentration, unlike fluctuations due to the intrinsic noise alone. Additionally, the mean concentration profile, averaged over both intrinsic and extrinsic noise, is elevated compared with the corresponding profile from a uniform sink distribution, showing that the classical homogenization approximation can be a biased estimator of the true mean.


I. INTRODUCTION
Transport processes in natural environments can involve an interplay between fine-scale disorder in the spatial domain within which transport takes place and randomness in the transport process itself. Theoretical models that seek to characterise outcomes in terms of means and covariances must therefore account for averages over the noise that is intrinsic to the transport process, and averages over the ensemble of random domains. Spatial averaging (via asymptotic homogenization or coarsegraining approximations) can be successful in capturing mean behaviour [1,2], but standard techniques often fail to quantify higher-order uncertainties. Here we use a simple reactive-transport problem to explore the relationships between intrinsic and spatial averages, and we present a hybrid homogenization method that predicts mean quantities and leading-order fluctuations due to the quenched disorder.
While interactions between intrinsic and extrinsic noise appear in applications ranging from gene expression to epidemic modelling [3][4][5][6], the problem we address is loosely motivated by physiology, an area in which predictive models are increasingly taking account of variability between (and within) individuals in order to inform personalized medicine [7]. In the placenta, maternal blood flows in a porous medium formed by a dense network of branches of villous trees, within which are capillaries containing fetal blood. Gas and nutrient exchange between mother and fetus takes place across the syncytiotrophoblast layer coating villous trees. Oxygen transfer between mother and fetus has previously been approximated using a simple one-dimensional model in which a chemical species moves via advection and diffusion past a spatially disordered array of point sinks [8,9], which take up the species via zeroth-order kinetics. The concentration of the substance post disorder average can (in general) be described using a homogenization approximation; fluctuations around the typical mean behaviour show long-range spatial correlation and have a structure and magnitude that is sensitive to both the statistics of the sink distributions and model parameters [8,9]. In some instances however, the fluctuations can become as great as the mean field itself and the homogenization approximation fails.
The present problem extends this work in significant respects. First, we treat the transport as a stochastic process, which enables us to exploit results derived for zero range processes [10][11][12][13]. Second, we assume the sinks operate via first-order kinetics and have variable strength rather than position. These features enable us to derive a hierarchy of descriptions that exploit the problem's multiscale structure, while remaining within a linear framework. Third, when the variance in sink strength is sufficiently small, we show how fluctuations due to the quenched disorder can be described analytically across a broad range of parameter space of our model (wider than that accessible to the direct method in [8,9]). These results can be used to examine systematic differences between averages over the sink strengths and averages over the intrinsic noise. These observations also illustrate differences between population-averaged results and outcomes predicted for an individual, and enable us to quantify the variability induced by the two distinct sources of disorder in the system.

II. MODEL
We frame our model in a generic manner in order to encompass both discrete and continuous transport processes. At the discrete level the model provides a simplified representation of (for example) the Brownian motion of a virus particle in a mucus film, with diffusive transport interrupted by adsorption at discrete sites on entangled macromolecules. At the continuum level, the model describes elements of the transport of a solute in a flow past an array of sinks, capturing some features of the porous medium encountered by maternal blood in the placenta, or airflow in a pulmonary acinus. Our main focus is on determining spatial characteristics of stationary-state particle distributions.

A. Model definitions and master equation
We consider M discrete sites, labelled i = 1, . . . , M , equally spaced in a domain of length L; see Fig. 1 for an illustration. The model describes one species of discrete particles moving in this domain. We write n i (t) for the number of particles located at the ith site at time t.
There is no upper limit on the number of particles that can reside at any site at any one time. The configuration of the system is determined by the site occupancies, written as n(t) = (n 1 , n 2 , . . . , n M ).
The model operates in continuous time. We assume there is an inflow of particles at the left boundary with constant rate α. Particles do not interact, so the influx is independent of the occupancy in the first site. In the bulk, each particle may hop one site to the right or left with rates p and q respectively. The total hopping rate from site i to i + 1 is then pn i , and that from i to i − 1 is qn i . Again there is no interaction between particles. Particles hopping to the right from the last site leave the system; the resulting outflow at the end of the chain is pn M .
Particles may also leave the system through a removal process at a subset of N sites that we call sinks; these are located at sites i 0 + ∆, i 0 + 2∆, . . . , i 0 + N ∆, where N ∆ + i 0 ≤ M . The integer ∆ is the sink-to-sink distance in units of sites. The particle removal rate at the j-th sink is S j n j∆+i0 , j = 1, . . . , N , if there are n j∆+i0 particles at the location of the sink. Using e i to denote the unit M -tuple with components e ij = δ ij , the transition rates in the model are therefore length of domain L number of sites M number of sinks N hopping rates p, q injection rate α mean uptake rate S0 physical distance between sites d = L/(M − 1) physical distance between sinks = ∆d number of sites/number of sinks The sink strengths S = (S 1 , . . . , S N ) will be treated as quenched random variables. They are independently drawn at the beginning, from a distribution f (S i ) with mean S 0 , variance S 2 0 σ 2 , and then remain fixed during the transport process.
We denote the conditional probability of finding the system in configuration n at time t, given a particular sink strength configuration S, by P (n, t|S). The particles hop according to a continuous-time Markov process with exponentially distributed waiting times between events. The time-evolution of the probabilities P (n, t|S) is governed by the master equation, with the transition rates as in (1). Using the model parameters p and q, and the inter-site distance d = L/(M − 1) and an inter-sink distance = d∆, we can identify a mean advection speed and diffusion coefficient as For later reference, we introduce a number of dimensionless parameters listed in Table 1. These include a Péclet number, based on the inter-sink distance, which characterises the relative strength of advection to diffusion, and a Damköhler number which characterises the relative strength of uptake to diffusion: 1. Illustration of the stochastic particle hopping model: qni, pni are the rates of hopping left, hopping right from site i; α is the rate of inflow at the left boundary; pnM is the outflow rate at the right boundary; Sjni is the removal rate at sink j (site i = j∆ + i0); and ni is the number of particles at the i-th site. ∆ is the number of regular sites between each pair of sink sites; in this figure ∆ = 3 and i0 = −2. The long-range dimensionless coordinate X ∈ [0, 1] spans the physical length L of the domain.
For the mathematical analysis in Sec. IV below we assume that the system contains a large number of sites and sinks (M, N 1). For later purposes, it is useful to introduce the inverse number of sinks, ε = 1/(N +1) 1. Our analysis applies for cases in which the sinks are sparsely distributed relative to the sites; we also introduce the ratio δ = 1/∆ ≈ N/M 1. We will refer to the noise due to the stochastic hopping as the intrinsic noise, and the disorder arising from the quenched sink strengths as the extrinsic noise. We write averages over the intrinsic noise (i.e., realisations of the stochastic hopping) as · · · I and averages over the extrinsic noise (i.e., the sink strengths) as · · · E .
For a fixed realisation S of the sink strengths we write ρ(t|S) ≡ n(t)|S I = n n(t)P (n, t|S).
This describes the (intrinsic) mean number of particles at the different sites at time t for fixed sinks S. Similarly we introduce an (intrinsic) covariance between the occupancies n i and n j , again for fixed sink strengths S, We write σ(t|S) for the resulting covariance matrix. The mean occupancies post intrinsic average in (5) can further be averaged over the extrinsic uncertainty. We use the following notation writing F (S) ≡ N i=1 f (S i ) for simplicity. The shorthand ρ(t) ≡ ρ(t|S) E is introduced for later convenience; overbars will be used to indicate averages over the extrinsic noise. The total expectation in (7) is an average over both sources of noise. Analogously, we can introduce and additionally the extrinsic covariance, The total covariance of n i (t) and n j (t) is then defined as where · · · I,E stands for the combined average · · · I E . After a modest amount of algebra one finds an expression of the law of total covariance. The first term in (10) is an average of the intrinsic covariance (6) over realisations of the sink strengths. The second term accounts for correlations between ρ i (t|S) and ρ j (t|S). These quantities are each obtained from averaging over the intrinsic noise only, but for a fixed realisation of the sink strengths. They will each depend on the sink strengths drawn, and can fluctuate together across realisations of S. Finally, we denote quantities in the stationary state of the dynamics by a superscript 'st'. For example, the stationary occupancies, averaged over the intrinsic noise, will be written as ρ st i (S). We will write ρ st (S) for the vector (ρ st 1 (S), . . . , ρ st M (S)).

III. NUMERICAL SIMULATIONS
In order develop a feeling for the behaviour of the model we first present numerical simulations. These are carried out in continuous time using the Gillespie algorithm [14,15]. We discuss two sets of simulations. The first set describes a case of densely spaced sinks, and is for a system of M = 10 sites with a sink at each site (∆ = 1, i 0 = 0). In the second set, sinks are more sparsely placed, specifically we use M = 100 sites, with sinks at every tenth site (∆ = 10, i 0 = 0). The remaining model parameters are S 0 = 1, p = 1, q = 0.5 and α = 100 in both cases.

A. Densely distributed sinks
We first consider a system with M = 10 sites, with a sink of strength S i = 1 at each site, resulting in Pe = Da = 2 3 for the above choices of p and q. There is no extrinsic disorder in this example. Removal is sufficiently rapid to prevent most particles from reaching ejection at the last site. Figure 2(a) illustrates the intrinsic stochasticity of the dynamics. We show a single realisation of the site occupancies n i (t|S) (solid lines), superimposed onto the mean occupancies ρ i (t|S), i = 1, . . . , 5 , obtained as an average of 10 4 independent runs. The intrinsic covariance matrix in the stationary state σ st ij is diagonal, see Fig. 2(b). We show in Section IV A that the occupancies n st i (S) and n st j (S) for i = j are independent random variables across realisations of the intrinsic noise whenever S is fixed.
In contrast, the total covariance in the stationary state will contain off-diagonal contributions when there is extrinsic uncertainty, as illustrated in Fig. 2(c) for a normal distribution of sink strengths S i with unit mean and variance 1/16. We note that a small proportion of the S i can be expected to be negative in this case; this does not have a significant bearing on the results in this example. The off-diagonal covariances imply spatial correlation between the intrinsic means of the occupancies at different sites across realisations of the quenched disorder.

B. Sparsely distributed sinks
A sparse distribution of sinks introduces a second length scale into the problem. This can be seen in Fig. 3(a), which compares the stationary mean occupancies for fixed sink strengths S i ≡ 1 (i.e., no extrinsic disorder) and normally distributed sinks (S i ∼ N (1, 1/16)). The parameters we use in this example result in Pe = 20/3 and Da = 400/3. The Damköhler and Péclet numbers are larger than in the previous example, i.e. sinkto-sink diffusion is weaker than before. Rapid removal at sinks again prevents most particles from crossing the whole domain, but the biased hopping is noticeable between each sink, with pronounced inter-sink staircases superimposed on a decaying profile of particle density. The total mean occupancy is slightly higher in the case of disordered sinks than in the case of constant sink strength S i ≡ 1, even though the number of sinks and their mean strength is the same in both examples; we explore the origin of this difference below. The intrinsic covariance in the case without extrinsic disorder (S i = 1 for all i) is again diagonal, see Fig. 3(b), whereas the total covariance with disordered sinks in Fig. 3(c) shows long-range spatial correlations and a multi-scale structure. The in- (c) Total stationary covariance σ tot,st (normalised as in (b)) for Gaussian sink strengths (unit mean, variance 1/16). Remaining parameters are M = 10, p = 1, q = 0.5, α = 100. Insets show the intrinsic variance σii (panel (b)) and the total variance σ tot ii (panel (c)) at each site. In order to show their raw magnitude, these are not normalised to unity. trinsic variance σ ii at the different sites shares the staircase structure of the mean occupancies, see the inset of Fig. 3(b). The total variance at the different sites has a striking non-monotonic form, as shown in the inset of Fig. 3(c). This indicates particularly strong variability immediately downstream of the first sink.
We now explore the origin of the long-range correla-

IV. ANALYSIS
We now proceed with a mathematical analysis of the model. An outline of our approach is illustrated in Fig. 4. We first briefly comment on the properties of the stationary distribution of the system (Sec. IV A). For a fixed realisation of the sink strengths we carry out an average over the intrinsic stochasticity and obtain the standard rate equations for the first and second moments of site occupancies; see Sec. IV B. These are ordinary differential equations (ODE), see also Fig. 4. Given that there are no interactions between particles (i.e., reaction rates are linear in the particle numbers), these equations close and do not involve higher-order moments. In a second step (Sec. IV C), and assuming a sufficiently large injection rate to ensure large particle occupancy at individual sites and a sparse sink distribution (δ 1), we take a continuum limit to derive a partial differential equation (PDE) for the mean occupancy, again for fixed realisations of the sink strengths. The PDE provides a continuum description of particle transport but retains a discrete representation of uptake at sinks. Then, assuming a large number of sinks across the domain (ε 1), we use a stochastic homogenization approach in Sec. IV D to obtain approximations for the total mean and covariance across the spatial domain. Whereas classical homogenization involves spatial averaging over a periodic microstructure to derive slow variation over macroscopic lengthscales, its stochastic analogue goes further by averaging over a disordered microstructure. In the present case, by assuming the disorder is weak, we will use the classical formulation as the starting point of a perturbation expansion in the small sink variance σ 2 . We validate these theoretical predictions against Monte Carlo simulations in Sections IV E-IV F. The range of validity of each of these approximations is assessed as a function of the input parameters of the model in Sec. IV G.

A. Stationary distribution, fixed sinks
The stochastic model, defined by the transition rates (1), is a variant of the open-boundary zero-range process (ZRP) [10][11][12][13]. It describes non-interacting particles, and includes particle removal dynamics. The stationary distribution of the open-boundary ZRP is a product distribution [10,11], i.e., in the stationary state the site occupancy numbers n st i , n st j are pairwise independent, and therefore uncorrelated. This distribution is independent of the initial condition, due to the ergodicity of the stochastic system. Following Levine et al.'s arguments [10,11], it can be shown these properties are left unchanged by the addition of particle removal through first-order sinks.  Using the results of [10], the stationary distribution of the model can be written in the form where the single-site marginal distributions are Poissonian. Their only parameters are the stationary mean occupancies ρ st i (S) = n st i |S I , for i = 1, . . . , M . We have Equation (11) can be evaluated if the stationary mean occupancies ρ st i (S) are known. Given the Poissonian nature of these distributions, we immediately conclude that the (intrinsic) variance at each site, for a fixed sample of the quenched disorder, equals the mean, σ st ii (S) = ρ st i (S). Furthermore, again for a fixed sample of the disorder, independence in the stationary state implies that the second moments factorize, n st i n st j I = n st i I n st j I , as earlier seen for example in Fig. 3(b). The total covariance in (10) finally becomes B. Exact equations for moments, fixed sinks The time-evolution of the means and covariances of the site occupancies n i can be derived directly from the master equation (2), see for example [16,17]. It is useful to define the M × M matrices A(S) and B(n, S) as for i, j = 1, . . . , M . We note that N k=1 δ i,k∆+i0 S k is the strength of the sink at site i, if there is one; this expression takes the value zero in absence of a sink at i. We also introduce the vector v with entries v i = αδ i,1 . Multiplying the expressions in (2) by n and summing over all configurations n yields (see Appendix A for details). Similarly, for a fixed sample of the quenched disorder the intrinsic covariances between the occupancies n i and n j satisfy [16,17]  When ρ st satisfies (15) in the stationary state, it is easily demonstrated that σ ij = ρ st i δ i,j satisfies (16). This is a consequence of Poissonian product form of the stationary distribution in (11).

C. Equations for moments in the continuum limit
We now consider the sites arrayed over a continuous spatial domain, and use (15) to derive a PDE for the first moment of the stochastic transport process at a fixed realisation of sinks. We approximate ρ i (t|S) by a continuous function C(x, t|S), where x ∈ [0, L] measures distance along the line of sites. One then has ρ i (t|S) = C(x i , t|S) for x i = (i − 1)d, i = 1, . . . , M . We retain the discrete locations of the sinks and introduce takes the form where the subscript t denotes a partial derivative and where we have used the definition (14) of the matrix A(S). We introduce nondimensional variables, denoted by asterisks, as where = d∆ is the physical distance between successive sinks and t 0 = 2 /D is the time scale of diffusion between sinks. The quantity C 0 drops out in (17), but it will be defined below. We also have where the factor is included to ensure that We substitute (18) and (19) into (17) and expand in δ ≡ 1/∆ = d/ 1. At a fixed number N of sinks, this is valid for large numbers of sites, M . We find This advection-diffusion-reaction equation is parameterised by Péclet and Damköhler numbers, defined in (4).
With multiple dimensionless parameters in the problem (Table I), it is important to distinguish carefully how each behaves when we take the limits of large site and sink numbers, while preserving low sink density. We analyse this a posteriori in Sec. IV G below. The equations at the inflow and outflow boundary sites differ from the bulk and must be treated separately. Under the scalings (18) the inflow boundary equation (15) becomes We can rearrange (3) to write the rate constants p and q as D/d 2 ± 1 2 (u/d) = (D/d 2 )(1 ± 1 2 Peδ) respectively. Expanding (21) in powers of δ and rearranging gives where we have introduced the concentration scale The time derivative is among the O(δ) terms in (22) that are neglected in the limit δ → 0; this implies that this approximation may not capture rapid variations in the inlet concentration at very early times. The leading-order inflow condition is obtained as Similarly, at the outflow boundary we take the final equation in (15), write it in terms of the nondimensional continuous variables, and consider only leadingorder terms in δ. We find The PDE system (20,24,25) provides a convenient route for approximating conditional means ρ st (S) and, from (13), the total covariance. Intersink transport is governed by the advection-diffusion equation (20); the inlet and outlet conditions are quasi-steady, with advection and diffusion contributing to the imposed flux ε in (24) and advection being sufficiently strong to enforce zero concentration at the outlet, see (25).
Since Pe and Da were defined with respect to the intersink distance in (4), they appear naturally as parameters in (20). The parameter ε appears in the domain length (0 ≤ x * ≤ ε −1 ) and the inlet flux. When Pe = Da = 0, the problem has steady diffusion-dominated solution C * = 1−εx * for which C * varies by O(1) across the whole domain, reflecting the balance between inflow and diffusion across all the sites implicit in (23). If we now assume ε 1 and consider increasing Pe and Da from zero, uptake first becomes important for Da = O(ε 2 ), when C * x * x * balances DaC * over a distance ε −1 ; advection first becomes important for Pe = O(ε), when C * x * x * balances PeC * x * over a distance ε −1 . It what follows we therefore formally consider the distinguished limit ε → 0 with Pe/ε and Da/ε 2 remaining O(1) ; these latter quantities are the Péclet and Damköhler numbers defined relative to the domain length L. This ensures that advection, uptake and diffusion are all of comparable magnitude.

D. Averaging over extrinsic noise
We now adopt a homogenization approach, spatially "smearing" the discrete sink locations and averaging over the sink strengths in (20,24,25). We write the sink strengths as S * i = 1 + σŜ i where theŜ i are independent random variables with unit variance. When σ is sufficiently small we may work withŜ i ∼ N (0, 1): a small number of sink strengths will then be negative, but this is not excluded by our formalism, and does not change the outcome; alternatively, for larger values of σ, we adopt a log-normal distribution.
In the stationary state and dropping asterisks from now on we must solve Splitting the concentration into its deterministic and fluctuating parts, Averaging (27) over the quenched disorder and using the fact that Ĉ while the residualĈ satisfieŝ When σ 1, we may obtain a leading-order approximation to C by neglecting σ 2 ĈŜ j E in (28), namely subject to (28b). We can use this to findĈ in (29), neglecting the O(σ) correction in that equation. We will then return to (28) to compute the O(σ 2 ) correction to C. The leading-order approximation for C in (30) contains a periodic array of sinks of fixed strength. At this level we have discarded the quenched disorder entirely. A classical two-scale asymptotic homogenization approximation may be adopted for this reduced problem [1]. The solution is represented as a series C = C (0) (x, X) + εC (1) (x, X) + ε 2 C (2) (x, X) + . . . , where we recall that ε = 1/(N +1) is the inverse number of sinks in the system. The short-range variable x is treated independently of the long-range variable X = εx. We recall that we have dropped asterisks before (26), and that x takes values in the interval [0, ε −1 ]; the variable X takes values in [0, 1]. A classical argument, described for example in [8], shows that the leading-order approximation depends only on X and satisfies X | x=0 = ε, C (0) | X=1 = 0. (31b) These are derived formally assuming Pe = O(ε) and Da = O(ε 2 ), which ensures a leading-order balance of advection, diffusion and uptake [8]. This linear problem can be solved directly, and has solution where φ ≡ Da + Pe 2 /4. The function C (0) varies smoothly over the length of the domain and provides a leading-order approximation to C in the limit of infinitely many sinks, ε → 0; higher-order terms C (1) , C (2) , . . . retain a dependence on x and capture the jump in the derivative of C across each sink.
Comparing (30) and (31a) illustrates the nature of the homogenization approach: the discrete sum Da N j=1 C(x)δ(x − j) has effectively been replaced by the continuous function DaC(x) in order to obtain the leading-order homogenized solution C (0) . This reflects the "smearing out" of the sinks, and captures the net effect of multiple sinks over long length scales. While this ansatz is appropriate for slowly-varying functions subject to periodic forcing, it cannot necessarily be adopted more generally. Fig. 5 illustrates, for four sets of (Pe, Da), how C (0) captures the sample mean over realisations of (26). The panels illustrate cases in which (a) strong uptake leads to rapid decay of the concentration field, (b) elevated advection displaces the concentration field towards the downstream end of the domain, (c) advection, diffusion and uptake are in balance across the domain, and (d) advection is dominant except in a narrow diffusive boundary layer upstream of the outlet. In panels (a,c), for which Pe 1, diffusion dominates at the inter-sink scale leading to smooth sample means. In contrast, when advection becomes significant at the inter-sink scale (as in Fig. 3(a), for which ε = 0.1), C (0) captures the solution averaged over sinks (with error of O(ε)) but fails to capture its internal staircase structure. Nevertheless, Figure 5 illustrates how (32), derived for Pe ∼ O(ε) and Da ∼ O(ε 2 ), provides a useful approximation across a wide range of nearby parameter space.

E. Quantifying extrinsic fluctuations
We now seekĈ. To solve (29), we neglect the O(σ) correction that is quadratic in the fluctuations and apply the homogenization ansatz to the term Da N j=1Ĉ δ(x − j), replacing it with DaĈ(x). The perturbations to sink strengthsŜ j vary abruptly from sink to sink so we retain their discrete form, using C ≈ C (0) to estimate the strength of each term. This yields the approximate sys-temĈ in 0 ≤ x ≤ ε −1 , subject to (29b). It is evident thatĈ involves multiple independent components, each forced by an individual sink. This formulation is related to the so-called Duhamel expansion in stochastic homogenization, for which formal convergence results are available [18]; similar approaches have been adopted in hydrology [19]. The Green's functionĜ(x, y) of (33, 29b) satisfieŝ and takes the form where We have introduced g(x) ≡ Pe sinh(φx) + 2φ cosh(φx). We writeĈ in terms ofĜ and form sums of independent random variables: where the integer i is such that i < x ≤ i + 1. The resulting sum depends on the slow variable X through the slowly varying functions C (0) and G ± . Combining the N independent random variables and approximating sums with integrals we obtain the approximate distribution ofĈ, in terms of the long-range coordinate X, aŝ Using (32, 36) and numerically integrating for different Pe and Da yields the variance predictions in Fig. 6. These show good agreement with Monte-Carlo estimates. When advection is strong, the variance increases with distance before falling to zero at the outlet. We can also use the approximation forĈ to compute the transverse covariances Cov T E (Ĉ(X)) ≡ Cov E (Ĉ(X),Ĉ(1 − X)) (derived in Appendix B). Fig. 6 confirms that the present analysis captures predictions of Monte Carlo simulations. Once again the correlation between mean sink occupancies varies smoothly over the entire length of the domain, despite the fluctuations being driven over much shorter lengthscales.

F. Influence of fluctuations on mean occupancies
We now return to C, using (37) to evaluate ĈŜ j E in (28). Using the fact that Because C 0 andĜ are smoothly varying functions, it is legitimate to employ the homogenization ansatz in the final step of (39). Thus a refined approximation of C is given by a homogenized version of (28a) as subject to (28b,c). This linear equation can be split into two parts, C = C (0) + σ 2 C (F ) , where C (0) satisfies (31), and the correction due to fluctuations in the sinks satisfies UsingĜ to solve for C (F ) we obtain, in long-range coordinates, It is straightforward to demonstrate that C (F ) (X) is nonnegative. Since C (0) (X) ≥ 0, the conditionĜ(x, y) ≤ 0 (illustrated in Fig. 5) is sufficient for the integral over the product in (42) to be non-negative. In (36a), the exponential is always positive, and each hyperbolic function in G − (x, y) is non-negative for 0 ≤ x, y ≤ ε −1 except for sinh(φ(y − ε −1 )) ≤ 0. Therefore G − (x, y) ≤ 0. Also, the relation (36b) only involves swapping x and y and an exponential factor, so G + (x, y) ≤ 0. Hence C (F ) (X) ≥ 0.
The correction is illustrated using the example in Fig. 3. We use ε = 0.1, implying that only limited accuracy can be expected of the homogenization approximation, and Pe = O(1) implying that the staircase structure appears at higher order in ε. In this case C (0) captures the decay in the mean concentration with distance reasonably well, while C (F ) captures the correct sense and magnitude of the correction due to fluctuations in sink strength.
Finally, to test how well this approach works for larger sink variances, we present simulations with log-normally distributed sink strengths, ensuring that S i > 0. Figure 7 compares simulations with σ 2 = 1 against the theoretical predictions of the mean (32), its correction (42) and the covariance (B3). The small-σ predictions of mean and variance provide surprisingly good approximations of both quantities. We now seek to understand in more detail the range of validity of the approximation.

G. Size of fluctuations
It is instructive to consider the outcome of the model in various regions of the space spanned by the parameters Pe and Da. Regime We can analyse the magnitudes of the contributions to the total covariance (13) from the intrinsic and extrinsic noise for each parameter regime. To do so, we estimate the magnitudes of C (0) andĜ by considering the dominant terms in governing equations (20) and (34) in the different regimes, and then use the estimates σ 2 C (F ) ∼ Da 2 σ 2 C (0) G 2 ± /ε (from (42)) and σ 2 Cov E ∼ Da 2 σ 2 [C (0) ] 2 G 2 ± /ε (from (B3)). The homogenization approximation fails when σ 2 C (F ) becomes as large as C (0) , or equivalently when the extrinsic fluctuations (measured by the size of their standard deviation) become as large as the mean concentration. We note also that Cov E should be multiplied by C 2 0 and ρ st i by C 0 to transform back to dimensionful variables; see (23). As we are only interested in the relative magnitude of mean and (co)variance we simply divide the mean by C 0 in Table II, where we summarise our results, assuming σ is no greater than O(1). The following picture emerges.
[D] When diffusion is dominant over uptake and advection, the extrinsic noise is always small because Da 2 σ 2 ε 3 . The correction to the total mean due to extrinsic fluctuations can be neglected. The variance is dominated by the intrinsic noise provided Da 2 σ 2 C 0 ε 3 . Fluctuations due to intrinsic noise are small compared to the mean occupancy (i.e. ρ st i /C 0 C (0) ) provided C 0 1.
[U] When uptake dominates advection (taking place over a length scale x ∼ Da −1/2 ), the correction to the total mean due to extrinsic fluctuations be-  comes significant for Da > ∼ ε/σ 2 , implying a breakdown in the homogenization approximation; the example in Fig. 7(a,c) sits at this threshold. The intrinsic noise becomes as large as the mean (i.e.
There are therefore two independent thresholds at which the system becomes strongly disordered, with the size of the parameter ε 1/2 σC 0 relative to unity determining which one dominates.
[A] When advection dominates, C (0) andĜ exhibit boundary layers of length x ∼ 1/Pe. Extrinsic fluctuations become dominant for Daσ > ∼ ε 1/2 Pe (the example in Fig. 7(b,d) sits just below this threshold) and intrinsic noise becomes as large as the mean for Pe > ∼ εC 0 .
These thresholds are illustrated in Fig. 8. The conditions on C 0 (see (23)) for intrinsic noise to be small compared to the mean can be re-expressed in terms of the parameters of the discrete model as the three conditions applying in the diffusion-, uptakeand advection-dominated regimes respectively.

V. DISCUSSION
We have investigated a model transport problem that incorporates both intrinsic noise associated with the underlying stochastic hopping process, and extrinsic disorder arising from variability in sink strengths. The former generates independent fluctuations in site occupancies, represented by a diagonal covariance matrix, typical of a ZRP. The latter generates long-range perturbations that can be correlated across the entire domain. We examined the case in which multiple sinks are distributed sparsely across the domain, allowing continuum multiscale approximations to be adopted. While it is natural to predict mean site occupancies using the ensemble-averaged sink strength (represented by the leading-order homogenized solution C (0) ), we found this to be a biased estimator of the true ensemble mean. This is because a locally elevated [diminished] sink strength leads to global reduction [increase] in concentration, including at the sink itself. This in turn leads to a net reduction in the average local uptake rate CS (represented by ĈŜ j E ≤ 0 in (28)). The homogenized solution therefore overestimates the uptake rate when there is variability in sink strength, and therefore underestimates the mean site occupancy.
We used stochastic homogenization to derive explicit predictions of the fluctuations arising from weak sink disorder, and validated the predicted covariance against simulations. The transport process has three competing physical effects -diffusion, advection, and uptake -and a relatively complicated interplay between these effects is observed.
The convergence of the homogenization approximation to the ensemble mean is parameter-dependent, weakening with increasing mean sink strength; i.e. with increasing Da in Fig. 8. The condition for homogenization to fail, σ 2 > ∼ max(ε/Da, εPe 2 /Da 2 ), can be expressed in terms of the parameters of the discrete model as which shows how the effects of disorder become important when the number of sinks falls and their strength increases. We estimated the relative magnitudes of intrinsic and extrinsic noise, showing how the former becomes prevalent as the inlet flux α diminishes (see (43)). Our analysis indicates that the parameter ε 1/2 σC 0 ∼ ε 1/2 σαM/(p + q) must be small compared to unity for intrinsic noise to dominate extrinsic noise.
In the present study we have not sought to describe the case of strong quenched disorder, defined by (44) and indicated by the shaded region in Fig. 8. We anticipate that individual realisations will deviate significantly from the ensemble average, making the system non-self-averaging in this parameter regime. Techniques from condensedmatter physics, such as the coherent medium approach and related methods [20,21], could be useful for estimating mean transport properties. Likewise we have not addressed time-dependent variations in detail, for which anomalous transport effects can be anticipated; this has been illustrated for a related chemical transport problem in the weak disorder regime [22], and framed as a continuous-time random walk [23].
Returning to one of our motivating problems, for oxygen transport in a placental subunit (a placentone), the Péclet number has been estimated to be of order 10 3 to 10 4 [24]. Taking the domain length L to be comparable to the path length (O(1cm)) from a spiral artery to a draining decidual vein, this implies Pe/ε > ∼ 10 3 , a regime in which advection dominates at the microscale. The spatial disorder of villous branches within the placentone will contribute to fluctuations in the concentration field induced by variability in uptake strength (as modelled here). Intrinsic noise due to small particle numbers can be expected to be negligible; however the influence of fluctuations in the flow field induced by the irregular geometry may be significant [25] and will be addressed elsewhere. An alternative application for which intrinsic and extrinsic noise may be of comparable importance concerns the motion of inhaled nanoparticles (such as viruses or drugs) through the mucus lining of a lung airway [26,27]: here predominantly diffusive transport may be mediated by trapping of particles by large mucin molecules. While the present model describes a limited number of features of such applications, it provides a framework for describing the magnitude, structure and influence of fluctuations.
The problem we have addressed has a number of obvious extensions, including spatially correlated or more densely distributed sinks, random sink locations, nonlinear kinetics and nonlinear hopping rates, higher spatial dimensions, etc. These extensions can be adapted to study specific applications in natural systems involving transport in the presence of spatial disorder. Of particular significance in terms of predictive modelling is understanding the nature and magnitude of the bias in the homogenization prediction. The present approach is a weak disorder expansion (see (44)) that allows the physical system to be described as a Gaussian process with slowly varying mean and spatial covariance. While this approach has wide applicability as a method of uncer-tainty quantification, alternative approaches are needed to address the strong disorder case in which extrinsic fluctuations appear at leading order. Since Cov E Ŝ k ,Ŝ l = δ k,l , the covariance simplifies to Cov E Ĉ (x),Ĉ(y) = Da 2 min(i,j) k=1 C (0) (εk) 2 G + (x, k)G + (y, k) + i k=j+1 C (0) (εk) 2 G + (x, k)G − (y, k) where the piecewise nature ofĜ takes care of the different sums. Then by approximating the above sums with integrals to leading order, we have the following expression for the covariance in long-range coordinates: Recall thatĜ varies by O(1) with respect to the slow variable X when Pe = O(ε), Da = O(ε 2 ).