## Abstract

Three stochastic-based methods are proposed for solving unsteady scalar transport problems in bounded, single-phase domains. The first (Method I), a local solution appropriate to problems having Dirichlet conditions, adapts a well-known local stochastic solution of a backward Fokker–Planck equation to scalar transport. Method II, a local solution applicable to Dirichlet, Neumann and/or mixed initial boundary value problems (IBVPs), and representing a time-dependent extension of a recently reported heuristic steady solution, provides a straightforward addition to the limited collection of techniques available for Neumann and mixed problems. This approach is shown to be equivalent to a long-standing, rigorous low-order solution and, in addition, allows development of a probabilistic-based analytical solution to Neumann problems, stated in terms of an exit probability. Method III, a global solution, likewise suitable for Neumann and mixed IBVPs, follows by combined application of domain boundary Taylor expansions and Method I. This approach is shown to be computationally equivalent to a global version of Method II.

## 1. Introduction

Stochastic methods have been used to study scalar transport for quite some time; significant, ongoing research can be found, for example, in the areas of ground water (Chevalier & Banton 1999; McKinley 1999; Nguyen 1999; Trantham & Durnford 1999; Grzywinski & Sluzalec 2000; Berkowitz & Scher 2002), atmospheric boundary layer (Gusey 1998; Reynolds 1999; Kurbanmuradov & Sabelfeld 2000; Kljun 2002; Ditlevsen 2003; Franzese 2003; Hseih *et al*. 2003) and marine (Reynolds 2002) transport. More recently, stochastic methods have found application in turbulent momentum, heat and mass transfer (Dekker *et al*. 1995; Papavassiliou & Hanratty 1997; Franzese *et al*. 1999; Mitrovic & Papavassiliou 2003; Mito & Hanratty 2003), and, for example, in the study of microstructure evolution during solidification (Charbon & LeSar 1997).

The theoretical basis of much of this work is founded on an equivalence that exists between continuum models describing scalar transport and, on a more fundamental level, stochastic differential equations that model the evolution of individual quasi-physical entities perpetuating transport (which we will refer to, interchangeably, as particles, random processes or random walkers (RWs)). For example, in the case where continuum-scale transport is determined by combined advection and Fourier diffusion, individual particle trajectories, **Χ**(*θ*) considered in this paper as moving in backward time *θ*, can be modelled using a stochastic differential equation of the standard Ito form,(1.1)where **b** is the drift velocity; is related to the diffusion tensor; **a** via , d*θ* is a backward time increment (d*θ*=−d*t*); and d**w**=d**w**(*θ*) is a multidimensional Weiner process. Other representations for d**Χ** are possible depending on the nature of the stochastic process modelled (e.g. Gardiner 1983).

Given (1.1), and assuming that the continuum scalar field *η*=*η*(**x**, *θ*) is twice differentiable in **x** and once differentiable in *θ*, we can apply Ito's formula to *η*(**Χ**(*θ*),*θ*) take expectations over all particle trajectories launched, and obtain a generic local solution of the form(1.2)where refers to the expectation taken with respect to the solution point , and where, as detailed in §3, the continuum local, advection and diffusion terms, having the form , sum to zero. Since the first two terms on the right-hand side of (1.2) are zero, the local stochastic solution for *η* is determined by the average of all boundary values sampled, plus the average of all initial values sampled, taken over the set of particles launched from .

Importantly, the equivalence between the continuum model describing scalar transport and the stochastic differential equation governing the motion of individual particles, e.g. (1.1), provides a Monte Carlo-based recipe for constructing stochastically based solutions for : launch sufficient number of RWs from , while appropriately accommodating the boundary and the initial conditions extant on the problem.

This paper focuses on passive scalar transport within single-phase regions. Three methods are proposed for solving unsteady problems on bounded domains, where the first applies to Dirichlet problems and the latter two are appropriate to Neumann or mixed Dirichlet–Neumann problems. The formulation of each method emphasizes heuristic and algorithmic reasoning, both as a means of illustrating the physical basis of each approach, and as a means of establishing clear, intuitive foundations for future development of numerical solution techniques.

The first method, Method I, represents a straightforward adaptation of a well-known local solution (Friedman 1975) to a backward Fokker–Planck equation, subject to Dirichlet boundary conditions; surprisingly, it does not appear that this solution has found wide application in scalar transport. Indeed, Method I proves important in several respects. First, the method accommodates ill-posed transport problems in which boundary conditions are unknown over a portion of the domain boundary (§3*c*). Second, the method provides the basis for development of Method III, a stochastic whole-field solution to mixed initial boundary value problems (IBVPs). Third, the method can be adapted to tackle nonlinear transport (§§3*d* and 5*c*).

Several significant results follow from the development and application of Method II, a local solution appropriate to unsteady Neumann, Dirichlet and mixed IBVPs. First, the method, which extends the steady heuristic approach of Grigoriu & Samorodnitsky (2003) to unsteady problems, is shown to be equivalent to Milshtein's rigorous low-order solution (Milshtein 1997). Second, with respect to unsteady Neumann problems, Method II provides the basis for a probabilistic-based analytical solution (§4*b*). The probabilistic solution, applied to diffusion-dominated problems in two limits, recovers, in both cases, corresponding continuum analytical solutions (§§4*b*(i),(ii)). Third, the latter comparisons provide consistent estimates for the particle reflection distance, *ϵl*_{0} (Grigoriu & Samorodnitsky 2003), used at Neumann boundaries, showing that *ϵl*_{0} should be of the order of the incremental diffusion length, , where 2*α* is the diffusivity and Δ*θ* is the backward time-step. This result is in turn consistent with and completes Milshtein's prescription for the thickness of his near-boundary region (which, as discussed in §4*a*, accommodates particle reflections). Fourth, Method II allows a straightforward demonstration that Method I can accommodate combined Dirichlet and homogeneous Neumann conditions (§5*a*), a feature that is crucial to the development of Method III.

Finally, Method III, which again applies to Neumann and mixed IBVPs, uses surface Taylor expansions to generate a system of equations involving unknown surface fluxes and unknown surface scalar magnitudes (extant on Dirichlet and Neumann boundaries, respectively), along with an associated set of subsurface scalar magnitudes. Application of a modified form of Method I, which again accommodates combined Dirichlet and homogeneous Neumann conditions, and which is used at each subsurface Taylor expansion location, then allows determination of subsurface scalar magnitudes, closing the problem. It is shown that by choosing the Taylor expansion distance, to be of the order of , Method III and a global version of Method II are computationally equivalent.

Given that stochastic methods offer a potentially powerful, mesh-free approach for studying scalar transport, particularly within small, spatially limited regions, it is clear that valid solutions in these cases require formulation of reliable boundary conditions. Thus, §6 of the paper considers an extended example of setting appropriate boundary conditions on a spatially focused stochastic solution. The example highlights the importance of properly determining a problem's structure prior to constructing a stochastic solution.

## 2. Preliminaries

### (a) Solution domains and boundary conditions

The methods to be developed apply to initial value transport problems on bounded domains. In all the cases, solutions are sought over a backward time-interval, , where backward time, *θ*, is related to forward time, *t*, via *θ*=*T*−*t*, , and where is a specified stopping time. We use backward time since it allows straightforward accommodation of non-homogeneous initial conditions, viz. these are effectively treated as Dirichlet conditions. Thus, known initial values are imposed on the final space-time slice, , where *D* is the spatial domain. Without loss of generality, all the three methods focus on two-dimensional transport within a rectangular domain *D*. The space-time solution domain, within which RWs are used to construct solutions, in all cases thus corresponds to a cylinder, denoted as *Q*, where (figure 1).

In Method I, Dirichlet conditions, or Dirichlet conditions combined with homogeneous Neumann conditions, must be applied on all the boundaries. The use of homogeneous Neumann conditions in the stochastic construction underlying Method I, while apparently not discussed in the literature, is detailed in §5*a*. Methods II and III allow for any combination of space and time-varying Dirichlet and Neumann conditions.

### (b) Random walks in fluid and solid phases

The space-backward time trajectory, , taken by any given RW is governed by an Ito stochastic differential equation of the form given by equation (1.1). In fluid-phase transport problems, the drift **b** assumes two distinct forms, depending on whether the flow is laminar or turbulent. In laminar flows, **b** corresponds to the reversed Eulerian velocity field, **b**=−**u**; by contrast, in turbulent transport problems, **b** is typically determined using Thomson's (Thomson 1987) well-mixed condition (e.g. Hseih *et al*. 2003), where again a sign change is required when backward solutions are sought. Considering the diffusion matrix, **a**, in laminar flows with isotropic diffusion, this term is simply related to the scalar diffusivity, 2*α*, via *a*_{ij}=2*αδ*_{ij}; while in turbulent flow, it is typically determined by the mean turbulent kinetic energy dissipation rate (Hseih *et al*. 2003).

For solid-phase transport, in cases where the solid undergoes time-dependent elastic or plastic deformation, the drift **b** again corresponds to the reversed Eulerian velocity field. Likewise, the form of the diffusion matrix, **a**, depends on the degree of anisotropy underlying diffusive transport; Carslaw & Jaeger (1959), for example, provide representative examples of anisotropic solid-phase diffusion.

Here, for simplicity, it is assumed that coupling between the scalar and the velocity fields is weak enough that the velocity field, **u**, can be computed independent of the scalar transport problem.

## 3. Transport problem subject to Dirichlet conditions: Method I

This section first reviews the standard local, backward-in-time stochastic solution for Dirichlet problems. Details of the stochastic construction and two simple examples are then used to illustrate the meaning of the solution. Finally, limitations associated with the solution are identified and remedies proposed.

The problem governing transport of scalar *η* is stated as follows:(3.1)(3.2)(3.3)where(3.4)and where *ϕ*(**x**) is the final condition on *η* and *g*(**x**, *θ*) is the time-varying Dirichlet condition. In fluid- or solid-phase scalar transport problems, *f*(**x**, *θ*) could represent a mass source or sink, for example, owing to chemical reaction, or a thermal source, for example, owing to viscous dissipation in fluids.

### (a) Local stochastic solution

The stochastic solution of (3.1)–(3.3), obtained at any time , , can be written in the form (Friedman 1975)(3.5)where estimates of expectations are given by(3.6)(3.7)and(3.8)In (3.6), *N* is the total number of RWs launched from , is the (random) exit time for the *i*th RW to exit through *δQ*, **Χ**(*τ*_{i})∈*δQ* is the corresponding exit point and *N*_{g} is the total number of trajectories exiting through *δQ*. Similarly, in (3.7), is the travel time for each RW reaching *δQ*_{f}, is the exit point and *N*_{ϕ} is the number of RWs exiting through *δQ*_{f}. Likewise, *τ*_{i} in (3.8) is the first exit time for trajectories passing out of either *δQ* or *δQ*_{f} (where in the latter case). A schematic of the solution construction is shown in figure 2.

### (b) Heuristic examples

As a means of gaining physical insight into the solution in (3.5)–(3.8), simple heuristic checks can be carried out. Thus, consider, for example, purely diffusive transport (**b**=**0**) within a slender domain *D*, in which the source term *f*=0 and in which the aspect ratio where *L*_{0} and *H*_{0} are the domain dimensions in the *x*_{1} and *x*_{2} directions, respectively. For times *θ* satisfying , and solution points lying in the range , most RWs will impinge on either lateral boundary, *x*_{2}=0 and *x*_{2}=*H*_{0} (rather than exiting through the final time slice, *D*×{*θ*=*T*} or through end boundaries, *x*_{1}=0 and *L*_{0}). Thus, *N*_{g}≫*N*_{ϕ}, so that *N*_{g}≈*N*. Considering the case where *g*(*x*_{1}, *x*_{2}=0, *θ*)=*g*(*x*_{2}=0)=*c*_{1} and *g*(*x*_{1}, *x*_{2}=*H*_{0}, *θ*)=*g*(*x*_{2}=*H*_{0})=*c*_{2} with *c*_{1} and *c*_{2} constant, transients die out for , and the leading-order transport equation assumes the form . The corresponding leading-order solution, having *O*((*H*_{0}/*L*_{0})^{2}) error, is *η*=(*c*_{2}−*c*_{1})*x*_{2}/*H*_{0}+*c*_{1}.

We verify that this solution is recovered via the stochastic formula in (3.5) by first noting that under the present set of conditions, the solution in (3.5) can be expressed in terms of exit probabilities, , where *ψ*_{1}(*x*_{2}) and *ψ*_{2}(*x*_{2}) are the probabilities of exit through *x*_{2}=0 and *x*_{2}=*H*_{0}, respectively, and where *ψ*_{1}(*x*_{2})+*ψ*_{2}(*x*_{2})=1. Here, the problem governing, for example, *ψ*_{1} is given by , with *ψ*_{1}(*x*_{2}=0)=1 and *ψ*_{1}(*x*_{2}=H_{0})=0 having the solution *ψ*_{1}(*x*_{2})=1−*x*_{2}/*H*_{0}. Thus, we find that the stochastic solution given by (3.5) is exactly equal to the continuum analytical solution.

A similar check of the limit where axial drift, given, for example, by , dominates diffusion, with associated boundary layer growth, , much smaller than *H*_{0}, shows that outside the boundary layers adjacent to *x*_{2}=0 and *x*_{2}=*H*_{0}, and for −*L*_{0}/*b*_{1,avg}, RWs released from any point (*x*_{1}, *x*_{2}) track to the inlet at *x*_{1}=0. In this case, again *N*_{g}≈*N*, and the time for particle exit through the inlet is given by *τ*≈*τ*(*x*_{1}, *x*_{2})=*x*_{1}/*b*_{1}(*x*_{2}); thus, . This, of course, corresponds to the continuum analytical solution obtained via the method of characteristics.

### (c) Limitations

#### (i) Unknown boundary conditions

This section and §3*d* briefly address two limitations associated with the Method I solution, namely the fact that many transport problems are subject to *a priori* unknown boundary conditions, and that transport problems are often nonlinear.

Considering the first limitation, we focus on the practically important example in which an outlet condition in a boundary layer development problem is unknown. The problem is shown schematically in figure 3, where the inlet and outlet correspond to the left and right boundaries, respectively. Under these circumstances, a certain fraction of the *N* RWs launched from a solution point will typically exit through the uncharacterized downstream boundary, .

This source of error can be circumvented by at least three approaches. The first, and often simplest, as illustrated in the extended example of §6, centres on establishing a realistic downstream exit condition, based on consideration of conditions external to the boundary layer. Alternatively, it is often possible to infer a homogeneous Neumann condition; in this instance, as detailed in §5*a*, the solution in (3.5) remains valid.

A second approach, which is the most rigorous, requires calculation of the exit probability, *ψ*(**x**, *θ*) for RWs leaving the cylinder *Q* through *δQ*_{e}. For example, considering the broad class of problems not subject to volumetric generation or consumption, we set *f*(**x**, *θ*)=0 in (3.1), and define the final and boundary conditions corresponding to (3.2) and (3.3) as follows:(3.9)(3.10)and(3.11)where, as shown in figure 3, *δD*_{p} is the union of the inlet, upper and lower boundaries of the solution domain *D*. By considering the stochastic solution given in (3.5)–(3.7), it is readily observed that this problem indeed provides the probability of exit through *δQ*_{e}; specifically, setting *ϕ=0* and *g*(**x**, *θ*)=0 on all boundaries save the outlet, where *g*=1 the sum in (3.7) is identically equal to the number of trajectory exits through *δQ*_{e}, divided by the total number of RWs launched. Thus, solving the relatively simple exit probability problem governed by (3.1) and (3.9)–(3.11), via, for example, analytical or numerical means, provides a detailed map of the conditional probability, , for trajectory exits through the outlet, as a function of the space-time solution point . In practice, one would probably define an *ad hoc* upper limit on *P*, above which stochastic solutions at any associated coordinate would not be accepted.

A third, heuristic, approach simply limits stochastic solutions to the region , where is the characteristic diffusion length and *L*_{0} is again the streamwise length of the solution domain *D*. Since the variance in the swarm of RWs launched from a solution point increases as 2*αθ*, the characteristic radius of the swarm at time *θ* is 2*αθ*^{1/2}. Thus, limiting solutions to the above streamwise range limits the number of trajectories exiting through *δQ*_{e}.

### (d) Nonlinear transport

A standard approach for treating nonlinear, time-parabolic transport problems is based on quasilinearization. The idea centres on approximating the *current* magnitude of any given transport parameter, in our case the diffusivity matrix, **a**(**x**, *t*)=**a**[*η*(**x**, *t*)], using the solution for *η*(**x**, *t*) obtained at the previous time-step,(3.12)In order to use quasilinearization, the transport problem must be parabolic in time, i.e. the problem must evolve from some initial state, *η*(**x**, *t*=0)=*ϕ*(**x**), **x**∈*D*.

Here, for a backward-in-time stochastic solution, the following sequential algorithm can be used. Starting from the initial forward time, *t*=0, discretize the forward time axis as *t*_{k}=*k*×Δ*t*, *k*=0, …, *M*, where Δ*t*=*T*/*M*, *T* is the desired total solution time and *M* is the number of time-steps. Given this discretization, and beginning with the initial *forward* time-interval (0, Δ*t*], apply the *backward* stochastic solution approach, Method I, in succession over the first (*k*=1), and then over each succeeding forward interval, [*t*_{k−1}, *t*_{k}). Thus, for example, the *final condition* to be sampled during the *k*th solution interval corresponds to *η*(**x**, *t*_{k−1}), as computed in the preceding time-step (where *η*(**x**, 0)=*ϕ*(**x**) for *k*=1). Likewise, only those boundary conditions extant on the space-time boundary increment Δ*Q*^{(k)}=*δD*×(*t*_{k}−1, *t*_{k}] are sampled over (*t*_{k}−1, *t*_{k}]. In addition, during each time-interval, (*t*_{k}−1, *t*_{k}], quasilinearization is used to evaluate all scalar-dependent properties. Alternative approaches designed to improve solution stability (at higher cost), for example, in which sampling occurs over, for example, *δD*×(0, *t*_{k}], are also possible.

## 4. Transport problems having mixed Dirichlet and Neumann conditions: Method II

This section and the next describe stochastic-based approaches for solving time-dependent transport problems subject to both Dirichlet and Neumann boundary conditions. Importantly, it appears that very few stochastic techniques have been developed for solving steady and time-dependent mixed boundary value problems (Milshtein 1997; Grigoriu and Samorodnitsky 2003).

The problem is depicted in figure 4 where, for illustrative purposes, Dirichlet conditions are imposed over *δD*_{1}, and Neumann conditions are imposed on *δD*_{2}. In general, *δD*_{1} and *δD*_{2} can each consist of finite sets of disconnected line segments (or in three-dimensional problems, disconnected surface areas). Finally, we only consider the case where diffusivity is isotropic, i.e. *a*_{ij}=2*αδ*_{ij}.

Method II extends Grigoriu & Samorodnitsky (2003) recently reported heuristic approach for solving *steady* mixed Dirichlet–Neumann problems. The idea is as follows: whenever an RW reaches a Neumann boundary, the RW is stopped, displaced normally into the solution domain *D* by a small fixed distance, *ϵl*_{0}, and then restarted. The Markov chain thus constructed is terminated when it finally reaches a Dirichlet boundary.

Here, we define a backward-in-time *n*-dimensional random walk, * Χ*(

*θ*) (where in the present two-dimensional case,

*n*=2) governed by the backward-in-time version of Ito's equation, (1.1). We again choose a backward-in-time approach since the method requires a known value for

*η*at some termination time. Here,

*η*is known on the final time slice,

*δQ*

_{f}=

*D*×{

*θ*=

*T*}, as well as on the Dirichlet portion, , of the boundary . The termination times in each case are thus and , respectively, where is a random variable corresponding to the time of first impingement on

*δQ*

_{1}and is again the chosen solution time.

We thus launch a stochastic process, **Χ**(*θ*), from solution point , lying within the cylinder *Q*, and imagine the process impinging on the Neumann boundary, *δQ*_{2}, at a series of space-time points , where, as a generalization of Grigoriu & Samorodnitsky (2003), we stop the process at each point of impingement and displace the process into *Q* according to(4.1)Here, the displacement is in the direction opposite the local outward normal, **n**(**Y**_{j}); *l*_{0} is a small constant distance; and *ϵ*≪1. The above relationship applies to the *j*th impingement, where 1≤*j*≤*q*, and where the RW is terminated when it either impinges on *δQ*_{1} or *δQ*_{f}; thus, *q* is the number of impingements on the Neumann boundary *δQ*_{2} that occur during the RW's journey to either *δQ*_{2} or *δQ*_{f}.

Given (4.1), we next write the following pairs of equalities, valid at each of *q* succeeding impingement points:(4.2)Following the *q*th impingement, the RW either reaches *δQ*_{1} or *δQ*_{f}, yielding, respectively,(4.3)or(4.4)Here, the first equality in each pair of (4.2) reflects application of Ito's formula to *η*(**x**, *θ*), followed by use of equation (3.1). The last term, *M*_{p}, in these equations, given by is a martingale; in taking the expectation over the swarm of RWs launched from , the sum of these terms, , is zero. The second, approximate equality in each pair of (4.2) represents a finite difference approximation to the surface flux, , at each impingement point.

Summing all the terms on the right- and left-hand sides of (4.2) and using the appropriate end value for *η*, given by (4.3) or (4.4), we obtain an important relationship between *η* at the solution point, , the set of fluxes sampled at the *q* impingement points on *δQ*_{2}, and the final known value of *η*,(4.5)where *β*(**Y**_{q+1}, *θ*_{q+1})=g(**Y**_{q+1}, *θ*_{q+1}) for final impingement on *δQ*_{1} and *β*(**Y**_{q+1}, *θ*_{q+1})=*ϕ*(**Y**_{q+1}) with , for final impingement on *δQ*_{f}.

The final local solution for is obtained by taking the expectation over a swarm of RWs launched from , yielding an expression that is a generalization of the Method I solution given by (3.5)–(3.8),(4.6)where *N* is the number of RWs launched, *N*^{(i)} is the number of boundary hits on *δQ*_{2} experienced by the *i*th RW and is the *p*th impingement point and time for the *i*th RW. In addition, (**Χ**(*τ*_{i}),*τ*_{i})∈*δQ*_{1} in the first term on the right-hand side of equation (4.6) is the first impingement point on *δQ*_{1} for the *i*th RW to reach *δQ*_{1}; *N*_{g} thus equals the number of RWs reaching *δQ*_{1}. Likewise, **Χ**(*τ*_{i})∈*δQ*_{f} in the second term is the impingement point on the final time slice, *δQ*_{f}, of the *i*th RW to reach *δQ*_{f}; *N*_{ϕ} is thus the total number of RWs to reach *δQ*_{f}. (Hence, *N*_{g}+*N*_{ϕ}=*N*). Finally, in the third term, and are the start and stop times corresponding to the (*p*−1)th and *p*th hits of the *i*th RW on *δQ*_{2}.

### (a) Validation of Method II

In order to validate the solution in (4.6), we compare it to a rigorous low-order solution reported by Milshtein (1997). However, we first note the consistency of Method II. Specifically, in the case where the entire boundary is subject to Dirichlet conditions, the last term on the right-hand side of (4.6) is zero, and we recover the solution obtained via Method I, (3.5)–(3.8).

Since Method II reduces to Method I when Dirichlet conditions are imposed, we focus on the case where only Neumann conditions are extant, with the source term, *f*(**x**, *θ*), zero, and the initial condition, *ϕ*(**x**), homogeneous. Thus, (4.6) simplifies to(4.7)

The time-dependent transport problem involving both Dirichlet and Neumann conditions can in principle be solved using *n*+2-coupled stochastic differential equations, where again *n* is the spatial dimension of the transport problem. Under the conditions defined above, this system assumes the form (Milshtein 1997)(4.8)where *I*_{A}(**Χ**) is the indicator function for the set *A I*_{A}(**Χ**)=1 when **Χ**∈*A*, *I*_{A}(**Χ**)=0 otherwise) and *μ*(*θ*) is the local time process associated with the stochastic process **Χ**(*θ*). In the general case (Milshtein 1997), the stochastic process *Y* accommodates scalar-dependent source terms within *Q* and/or on *δQ*. Likewise, process *Z* accommodates scalar-independent source terms in *Q* and/or on *δQ*. In all cases, the initial values of *Y* and *Z*, evaluated at the solution point , are specified as 1 and 0, respectively. See Gihman & Skorohod (1972) and Milshtein (1997) for further details and for the general form of the coupled stochastic differential equation (SDE) system as well as the associated general probabilistic solution.

Milshtein (1997) derived both high- and low-order approximate solutions to the general stochastic system, appropriate for numerical treatments. Thus, following Milshtein (1997), we define a near-boundary region, , consisting of all **x** lying within a distance, *C*_{1}*r*, of the domain boundary *δQ*, where *C*_{1}*r* is measured in the local normal direction of *δQ*. Here, *r* is small and *C*_{1} is a positive *O*(1) quantity. Given , define as the annular layer adjacent to (and including) the space-time boundary .

Milshtein proved that once the process **Χ** enters the near-boundary region, , to an order of accuracy, *O*(*r*), the following formulae can be used to estimate the non-zero stochastic differentials in (4.8):(4.9)where again the differential d*Y*=0, and *q* is an unspecified *O*(1) positive quantity. Note that **n** is here defined as an inward unit normal; thus, as dictated by Milshtein (1997), and consistent with Grigoriu's prescription (Grigoriu & Samorodnitsky 2003), the displacement, Δ**Χ**.**n**=*qr*, of the random process, **Χ**, is away from the boundary, *δD*, into *D*. Also note that the time-step is given by Δ*θ*=*r*^{2}, where a larger time-step is permitted outside of .

Referring to the expression for d*Z* in (4.8), and to the near-boundary approximation, Δ*Z*, in (4.9), it is clear that for any given RW, say the *i*th RW launched from a solution point , that integration yields(4.10)where the terms and again represent the *p*th impingement location and time, respectively, and where the second term on the left-hand side again, is specified as 0. Taking the expectation over *N* RWs launched from , and using the general stochastic solution given in Milshtein (1997), we finally obtain Milshtein's *O*(*r*) solution(4.11)Comparing Milshtein's solution, (4.11), with the solution from Method II, (4.7), it is clear that the latter is equivalent to the former. We note that derivation of Method II, as given above, uses a fixed time-step, Δ*θ*, whereas Milshtein's approach allows for differing time-steps inside and outside the near-boundary region. In the case where a fixed time-step is used in the latter approach, the following correspondences between parameters can be identified: Δ*θ*_{II}=Δ*θ*_{M}*=r*^{2} and *ϵl*_{0}=*qr*, where the subscripts ‘II’ and ‘M’ refer to Method II and ‘Milshtein’, respectively. Finally, as discussed by Milshtein (1997), an *O*(*r*^{2}) near-boundary scheme can be generated by defining a near-boundary region as , where the thickness *δ* is chosen to be proportional to *r*^{2}. In this case, *r* in all of the differentials given above is replaced by *δ*.

### (b) Derivation of probabilistic-based analytical solution

In this section, we use Method II to derive a probabilistic-based analytical solution applicable to time-dependent Neumann problems. We then show that the probabilistic solution matches corresponding continuum analytical solutions in two limits (§§4*b*(i),(ii)). The latter analyses also provide detailed expressions for the unspecified parameters, *ϵl*_{0} and *qr*, introduced in Grigoriu & Samorodnitsky (2003) and Milshtein (1997) formulations, respectively. With regard to numerical solutions of problems that are beyond probabilistic analytical solutions, both examples indicate consistent estimates for these parameters.

Importantly, the *approach* used to obtain the probabilistic analytical solution appears to be general enough to be applicable to a range of transport problems subject to Neumann conditions. The two examples considered in §§4*b*(i),(ii) clearly indicate the suitability of this approach to diffusion-dominated problems. In problems where drift is important (and boundary injection/suction does not occur), the key assumption underlying the derivation of the probabilistic solution (i.e. that RWs, upon reaching a Neumann boundary, remain close to the boundary through the remainder of the solution) probably remains valid. Thus, for problems subject to advection, it likewise appears that the probabilistic-based solution still applies. Future work will consider this important question.

For simplicity, we assume diffusion-dominated transport and, in this section, limit solution points, (*x*_{1}, *x*_{2}) to the region and ; in this case, most RWs will not reach the lateral boundaries nor the upper boundary, so the problem can be taken as one dimensional (in the *x*_{2} direction) over a semi-infinite domain. In addition, we again assume a homogeneous initial condition and that sources, *f*(**x**, *θ*), are absent. Although the boundary flux, , is taken to be constant throughout, for linear problems, the stochastic solution obtained can be used to construct, via Duhamel's method (e.g. (Carslaw & Jaeger 1959)), solutions to problems in which the boundary flux varies with time. This is briefly discussed in §4*b*(iii).

To begin the derivation, discretize the *forward* time axis into *M* equal increments, with Δ*t*=*T*/*M*; in addition, for notational simplicity, let the backward solution time . Next, launch *N* RWs from the solution point . Over the *j*th forward time-interval, (*t*_{j−1}, *t*_{j}] (where *t*_{j}=*j*Δ*t*, *j*=1,2, …, *M*), let the total number of RWs that impinge on the boundary, *x*_{2}=0, be denoted as *n*_{j}. The stochastic solution in (4.7) can then be stated as(4.12)where is the constant imposed flux.

Next, we recognize that once an individual RW, travelling in backward time, has reached *x*_{2}=0, say over the *k*th backward time increment Δ*θ*_{k}=[*θ*_{k−1}, *θ*_{k}) (where *θ*_{k}=*k*Δ*θ*=*kT*/*M*, *k*=1,2, …, *M*) then owing to the small displacement *ϵl*_{0} away from *x*_{2}=0, the RW will, with probability approaching 1, impinge on the boundary again during the next and all subsequent backward time increments, Δ*θ*_{k+1}, Δ*θ*_{k+2}, …, Δ*θ*_{M}. Thus, defining as the number of RWs reaching the boundary *x*_{2}=0 for the first time (i.e. without any prior boundary hits) during the *k*th *forward* time-interval, Δ*t*_{k}, the following approximate equalities can be written as(4.13)

Thus, using (4.13) in (4.12) yields(4.14)

Next, note that (*M*−*k*)=(*T*−*k*Δ*θ*)/Δ*θ*, *k*=1,2, …, *M*, and recognize that the following relations can be written as:(4.15)where *ψ*_{0}(*x*_{2}, *θ*) is the probability that RWs launched from (*x*_{2}, *θ*)=0 exit through *x*_{2}=0 over [0, *θ*). Thus, using (4.15) in (4.14) gives(4.16)or(4.17)where *θ*=*T*.

The first bracketed term in (4.17) simplifies to *θ*[*ψ*_{0}(*x*_{2}, *θ*)−*ψ*_{0}(*x*_{2}, 0)]. Focusing on the second bracketed term in (4.17), we rewrite this asand use this in (4.17) to finally obtain(4.18)where a small-order term, *θψ*_{0},_{θ}|_{θ}Δ*θ*, has been neglected relative to the integral shown.

Importantly, and as noted, (4.18) comprises what appears to be a fairly general probabilistic-based analytical solution to diffusion-dominated problems subject to Neumann conditions. The following subsections compare this solution in two limits against two known analytical solutions.

#### (i) Comparison of probabilistic and continuum analytical solutions: X_{T}≫1

Both probabilistic- and continuum-based solutions are conveniently expressed in terms of , as shown in this section, in the limit *X*_{T}≫1, the probabilistic solution in (4.18) is asymptotically equal to the exact continuum solution.

For one-dimensional, diffusion-dominated transport over a semi-infinite domain, the exit probability, *ψ*_{0}, is governed by(4.19)with *ψ*_{0}(*x*_{2}=0, *θ*)=1, *ψ*_{0}(*x*_{2}→∞, *θ*)=0 and *ψ*_{0}(*x*_{2}, *θ*=*T*)=0. In order to obtain a solution, the problem is restated in terms of forward time using the substitution *θ*=*t*−*T*. The result, stated in terms of backward time, is given by(4.20)where .

Inserting (4.20) into (4.18), and focusing on the limit *X*_{T}≫1, we then obtain(4.21)where again, for notational convenience, we set . Note that the order of error in this solution is .

The continuum analytical solution to the forward-time, one-dimensional diffusion problem, evaluated at *t*=*T* (corresponding to *θ*=0 and subject to the conditions at *x*_{2}=0, *η*→0 as *x*_{2}→∞, and a homogeneous initial condition, is given by(4.22)

Comparing the asymptotic stochastic solution given by (4.21) with the analytical solution in (4.22), it is seen that by properly choosing the displacement, *ϵl*_{0}, the solutions match; thus, we choose(4.23)where is the characteristic diffusion distance associated with the time-step Δ*θ*.

Practically speaking, (4.23) is important since it provides an explicit relationship for estimating *ϵl*_{0}, or equivalently, *qr*. Specifically, since Method II is equivalent to Milshtein's *O*(*r*)=*O*(*ϵ*) scheme (Milshtein 1997), for numerical solutions based on Method II, the choice , for the displacement, where(4.24)and where *c*_{0}=*O*(1), ensures that the error associated with Neumann boundary–particle interactions remain *O*(*ϵ*). A similar relationship is obtained in the next example (§4*b*(ii)).

#### (ii) Lumped model:

Here, as a second application of the probabilistic solution in (4.18), we assume that the diffusive time-scale across domain *D*, , is short relative to the time, *θ*; that is, diffusion across the layer is rapid enough to allow a lumped continuum analysis. We again limit the problem to one-dimensional transport by assuming that *H*_{0}/*L*_{0}≪1, and likewise assume that over , a constant flux, , exists on both boundaries *x*_{2}=0 and *H*_{0}. The last assumption allows straightforward application of (4.18). Again, for linear problems subject to a time-varying flux, the last constraint is not restrictive since a stochastic solution constructed for the constant flux case can be incorporated into a Duhamel solution.

Define again the solution point and time as , and discretize the backward time axis into *M* equal increments, Δ*θ*. Once a RW is launched, its first impingement on one of the boundaries, *x*_{2}=0 or *x*_{2}=*H*, will occur in a time of the order of , where *τ*_{H}≪*T*. Since the same flux acts on both boundaries, the interaction of the set of *N* RWs with both boundaries (launched from ) is equivalent to interaction with a single boundary. In other words, for each time-interval, [*j*Δ*θ*,(*j*+1)Δ*θ*), *j*=0,1, …, *M*−1, all of the RWs reaching either *x*_{2}=0 or =*H*_{0} can be binned together; by this approach, the analysis leading to (4.18) still holds.

Thus, we replace *ψ*_{0} in (4.18) with *ψ*, the probability of exit through either *x*_{2}=0 or *H*_{0}, where *ψ* is governed by(4.25)with *ψ*(0, *θ*)=1, *ψ*(*H*_{0}, *θ*)=1 and *ψ*(*x*_{2}, *T*)=0, where *x*_{2}∈[0, *H*] and *θ*∈[0, *T*). Since , then, for , (4.25) simplifies to , with the corresponding solution, *ψ*=1. Using this solution in (4.18) finally yields(4.26)

The analytical solution follows at once via a lumped analysis(4.27)Again equating the stochastic and analytical solutions in (4.26) and (4.27), we obtain a relationship between *ϵl*_{0} and the time-step Δ*θ*,(4.28)

Consistent with the result in (4.23), the required displacement, *ϵl*_{0}, is proportional to , where the *x*_{2} length-scale, *x*_{s}, in the present case is *H*_{0}. Again, this result suggests that numerical treatments using , preserve *O*(*ϵ*) error for particle interactions with Neumann boundaries.

#### (iii) Time**-**varying Neumann conditions

In linear problems, the above results can be extended to allow stochastic solution of time-dependent Neumann problems. The approach uses a Duhamel integral solution, which, for example, in the problem discussed in §4*b*(i), assumes the form (Carslaw & Jaeger 1959)(4.29)where *ζ*(*x*_{2}, *t*) is the solution given by (4.21), stated in terms of forward time, *t*=*T*−*θ*, with having magnitude 1. Here, the time-varying boundary flux is given by .

## 5. Global solution of mixed Dirichlet–Neumann problems: Method III

This section proposes a second, *global* approach for tackling mixed initial boundary value scalar transport problems. For simplicity, we initially describe a whole-time-domain formulation and subsequently discuss a time-partitioned solution approach, suitable for both linear and nonlinear problems.

As discussed in §5*b*, the proposed approach is computationally equivalent to a global application of Method II. Conceptually, however, an essential difference, rooted in the treatment given random walk interactions with Neumann boundaries, differentiates the techniques. In Method III, non-homogeneous Neumann boundaries are made homogeneous and a slightly modified version of Method I, allowing standard reflections from these boundaries, is then applied (with Neumann conditions entering via surface Taylor expansions). In Method II, as described earlier, RWs sample non-homogeneous Neumann boundaries during the process of reflection.

As a preliminary, we state the stochastic solution of Method I in global form. Thus, discretize the spatial domain *D* into *N*_{0} sub-areas (sub-volumes in three-dimensional problems), the surrounding boundary *δD* into *N*_{s} discrete lengths (areas), and again, for simplicity, divide the time-interval [0, *T*) into *M* equal increments, Δ*θ*. Finally, define *N*_{M}=*N*_{0}.*M* and *N*_{T}=*N*_{s}.*M*, and assume that at each time-step, discrete values of both *η* and the source strength, *f*, are associated with each of the *N*_{0} sub-regions in *D*.

The local solution in (3.5) can then be stated in global form as(5.1)where the *N*_{M}×1 vector {**η**} is defined as {**η**}=[{**η**^{(1)}},{**η**^{(2)}}, …, {**η**^{(M)}}], and where each {**η**^{(j)}} represents the set of *N*_{0} scalar values extant over *D* at *θ*_{j}=(*j*−1)Δ*θ*. Likewise, {**ϕ**} is the discrete *N*_{0}×1 vector of initial values on *δQ*_{f}, {**g**} is the *N*_{T}×1 vector of boundary values on *δQ*, and {**f**} is the *N*_{M}×1 vector of source magnitudes over *Q*. Finally, the matrices **G**, **H** and **M**, having dimensions *N*_{M}×*N*_{T}, *N*_{M}×*N*_{0} and *N*_{M}×*N*_{M}, respectively, are determined by straightforward application of the local solution represented by (3.5), carried out at all *N*_{0} spatial evaluation points in *D*, during each instant, *θ*_{j}=(*j*−1)Δ*θ*, *j*=1,2, …, *M* spanning [0, *T*).

Shifting now to formulation of Method III, two steps underlie the approach. First, a vector of surface scalar magnitudes, {**η**_{s}}, extant on *δQ*=*δD*×(0, *T*], is related to associated sets of normal surface flux magnitudes, , and subsurface scalar magnitudes, , via Taylor expansions(5.2)Second, Method I is applied to generate a second linear system relating and {**η**_{s}},(5.3)where the overhead carets emphasize that these matrices are generated by launching RWs only from the set of *N*_{s}.*M* subsurface locations associated with {**η**}. (Note, the dimensions of each term in (5.2) and (5.3) are as follows: , , , , , , , and ).

Next, combining (5.2) and (5.3) and solving for {**η**_{s}} yields a closed system of the form(5.4)For example, in the case where the Neumann and Dirichlet portions of *δD*, *δD*_{2} and *δD*_{1}, respectively, are fixed in time, the vector {**η**_{s}} contains *N*_{1}.*M* known and *N*_{2}.*M* unknown elements, where *N*_{1} and *N*_{2} are the respective number of discrete elements comprising *δD*_{2} and *δD*_{1} (with *N*_{1}+*N*_{2}=*N*_{s}). Likewise, the vector contains *N*_{2}.*M* known and *N*_{1}.*M* unknown elements. Thus, (5.4) can in principle be solved for all (*N*_{1}+*N*_{2}).*M*=*N*_{T} unknowns. Given this solution, which in itself may be useful, for example, in inverse problems, a whole-field or local solution can then be obtained, for example, by application of Method I.

### (a) Homogeneous Neumann conditions and Method I

In order to generate, via Method I, the matrices , and in (5.3), homogeneous Neumann conditions must be imposed on *δQ*_{1}=*δQ*_{2}×[0, *T*), the Neumann portion of *δQ*. (No provision exists for treating non-homogeneous Neumann conditions via Method I.) Although apparently not discussed in the literature, the stochastic solution in (3.5) in fact allows imposition of zero-flux boundary conditions. This is seen most easily by considering Method II. Specifically, in the case where the boundary *D* is subject to both Dirichlet and zero-flux Neumann conditions, and where an individual RW impinges on the Neumann boundaries *q* times prior to reaching a Dirichlet boundary, the formula in (4.5) assumes the form(5.5)where again, *β*(**Y**_{q+1}, *θ*_{q+1}) is the value of *η* at the location, (**Y**_{q+1}, *θ*_{q+1}), where the RW first reaches either a Dirichlet boundary or the final time slice, and where all fluxes, , are zero. Thus, it is clear that imposition of zero-flux conditions on *δQ*_{2}=*δD*_{2}×[0, *T*) does not alter the Markov character of the individual random walks used to construct a local solution; on taking expectations at each of the space-time solution points associated with , one recovers (5.3) (or equivalently, considering individual elements of , equation (3.5)).

### (b) Equivalence of Methods II and III

Method III and a global version of Method II, both applied to generating the system in (5.3), are compared in table 1. As shown, the computational cost, indicated by consideration of the first three rows in the table, is comparable for both approaches. As also shown, respective errors incurred by RW interactions with Neumann boundaries are *O*(Δ*n*) and *O*(*ϵl*_{0}), where Δ*n*, the Taylor expansion distance, is a truncation error. Thus, by choosing Δ*n* to be of the order of *ϵl*_{0}, it is clear that, computationally, both approaches are equivalent.

### (c) Time-sequential algorithm for nonlinear and linear problems

A time-sequential algorithm, applicable to either nonlinear or linear problems, represents a combination of the nonlinear algorithm described in §3*d* and the solution described in §5. In particular, we use the same sequential approach described in §3*d*: Method I is applied, in turn, to successive increments, Δ*Q*_{j}=*D*×(*t*_{j}, *t*_{j}, +Δ*t*] of the space-forward time domain *Q*=*D*×(0, *T*]. Over each such increment, the method in §5*a* is then applied.

Thus, considering in detail the solution over Δ*Q*_{j}, we again have two independent relationships involving the vectors of instantaneous near-surface and surface scalar magnitudes, and , respectively, and the corresponding vector of surface fluxes, ,(5.6)and(5.7)where (5.6) represents the instantaneous system obtained via surface Taylor expansions and (5.7) is the instantaneous version of (5.3), the system obtained by executing the backward-in-time stochastic solution of Method I. Again, the *forward* time increment over which RWs sample boundary values of *η* is where random walk paths are determined using the backward-time approach of Method I. In addition, the *final* condition on each such solution corresponds to the solution obtained over the previous forward time-step; thus, as shown above in (5.7), {**ϕ**} in (5.3) is replaced by {**η**^{(j−1)}}. Finally, the matrices and are those obtained via the Method I stochastic solution over (*t*_{j}, *t*_{j}+Δ*t*]; in nonlinear problems, the transport properties are again determined by quasilinearization.

## 6. Boundary conditions and spatially focused solutions

Stochastic solution methods provide a potentially powerful tool for investigating scalar transport. The purpose of this section is to illustrate the importance of determining a problem's physical structure prior to setting up a stochastic solution. This initial assessment serves two important purposes. First, it allows proper formulation of an efficient, spatially focused stochastic solution. Second, it provides reliable boundary conditions on the stochastic solution.

For purposes of illustration and owing to its relatively rich structure, we focus on solution of the periodic Graetz problem in which the scalar, *η*, evolves within a fully developed laminar flow, subject to both a time-periodic variation at an inlet, *δD*_{i}, and position-dependent variations on upper and lower boundaries of a wide rectangular duct. The assumption of fully developed laminar flow is used for illustrative purposes only and is not limiting; as long as coupling between the scalar and the velocity fields is weak, the velocity field, **u**, whether developing or fully developed, laminar or turbulent, can be computed independent of the scalar transport problem.

For simplicity, a homogeneous initial condition is assumed and the volumetric source term, *f*(**x**, *θ*) in (3.1), is zero. Finally, non-dimensional terms are denoted with a tilde.

The non-dimensional form of the problem considered is stated as follows:(6.1)(6.2)(6.3)where(6.4)is the Peclet number based on the input disturbance wavelength *λ*_{0}=*U*_{0}/*ω*_{0}, *U*_{0} is the characteristic speed of the flow, *ω*_{0} is the inlet disturbance frequency, and where , and . Two other dimensional parameters arise, the streamwise length-scale, *L*_{0}, and the height, *H*_{0}, of the channel; throughout, we assume that *H*_{0}/*L*_{0}=*O*(1).

It turns out that the Peclet number, Pe_{λ}, features prominently in this problem. In particular, Pe_{λ} determines the structure of the core region outside the developing boundary layers, where Case I described below corresponds to the limit Pe_{λ}≫1, and where Case II applies when Pe_{λ}≪1. Importantly, in both instances, under the assumption that the boundary layer thickness remains significantly smaller than the cross-stream dimension, *H*_{0}, spatially focused stochastic solutions of the boundary layer regions can be set up. By contrast, when boundary layer growth is rapid enough such that *δ*_{max}/*H*_{0}=*O*(1), use of a focused stochastic solution domain is obviated. Here, in order to deal with error associated with a typically unknown downstream boundary condition, consideration of the exit distribution for Brownian trajectories from the solution domain becomes particularly important (see §3*c*(i)). In the following, we define , where *δ*(*L*_{0}) is the characteristic boundary layer thickness at *L*_{0}.

### (a) Case Ia: Pe_{λ}≫1; ;

When Pe_{λ}≫1, the structure of the problem is as depicted in figure 5. This subsection focuses on the non-boundary layer domain labelled region I, while §6*b* treats region III. In both cases, it is assumed that the boundary layer thickness at the end of the solution domain, *δ*_{0}=*δ*(*L*_{0}), has not grown to a significant extent relative to the cross-stream dimension, *H*_{0}. While one can clearly solve for the core scalar field via a non-stochastic numerical attack, in order to clearly expose the problem's structure, an analytical approach is preferred. Thus, as depicted in figure 5, depending on the length of the duct, *L*_{0}, one, two or three distinct regions within the core (non-boundary layer) region can be identified. In the first (region I), which begins at the point where boundary layer growth is initiated, *x*_{1}=0, the axial length-scale, *L*_{1}, corresponds to the wavelength, *λ*_{0}=*U*_{0}/*ω*_{0}, of the input disturbance in *η*; likewise, the inverse disturbance frequency, , determines the associated time-scale.

Thus, the leading-order forward-time equation, accurate to , governing the core flow in region I follows directly from (6.1)–(6.3),(6.5)The near-inlet scalar field thus propagates in wave-like fashion, with diffusive smearing representing a second-order effect. The leading-order solution, obtained via characteristics, is given by(6.6)

Importantly, for stochastic solutions limited to or encompassing 0<*x*_{1}*λ*_{0}, (6.6) provides an outer condition, accurate to , on the boundary layer adjacent to *x*_{2}=0 (and for that adjacent to *x*_{2}=*H*_{0}).

### (b) Case Ib: ; ; *x*_{1}/*λ*_{0}

The time-scale now corresponds to the characteristic diffusion time, , required to smooth peaks in the input variation in *η*, while the length-scale still equals *λ*_{0}. This case again describes the conditions extant in core region III of figure 5, where the left boundary, , corresponds to the product of *U*_{0} and *t*_{s}. Thus, rescaling (6.1) and (6.2), we obtain the equation governing in region III,(6.7)where the cross-stream diffusion term has been neglected. Thus, the core solution, which again would serve as an approximate outer boundary condition on a stochastic-based solution of the region III boundary layer development problem (accurate to , is(6.8)where is the ambient scalar magnitude within the flow.

Regarding the core problem in region II, treatment of this part of the problem requires inclusion of both the time and advective rates of change terms and the axial diffusion term. If a stochastic-based boundary layer solution is desired here, the simplest approach would use a numerical or analytical solution of the core, initiated from the inlet, followed by a focused stochastic solution over the streamwise range of interest.

Finally, in closing discussion of Case I, we note that the assumption that the boundary layer remains thin over the axial length *L*_{0}, is precisely stated asFor , this assumption is satisfied under a variety of conditions, for example, when *λ*_{0}/*H*_{0}≪1 or when *λ*_{0}/*H*_{0}≪*O*(1), (where, as stated above, *L*_{0}/H_{0}=*O*(1)).

### (c) Case II: ; ; *x*_{1}/*λ*_{0}=*O*(1)

This case applies in the near-inlet region and arises under conditions where the axial diffusion length-scale, *x*_{diff}=(2*α*/*ω*_{0})^{1/2}, is much larger than the disturbance wavelength, *λ*_{0}. In particular, taking the ratio of these scales shows that . Under such conditions, typified, for example, by high-frequency input conditions, the wavy near-inlet variation in *η* becomes diffusively smeared over a few wavelengths.

Recognizing that *λ*_{0} still represents the axial length-scale and that the appropriate time-scale is *t*_{diff}, we obtain the dimensionless equation for within the core(6.9)where again the cross-stream diffusion term is neglected. The leading-order solution, accurate to *O*(Pe_{λ}), is given by(6.10)where and . In analogy with Stoke's solution for viscous diffusion of vorticity near an oscillating flat wall (Panton 1996), the present solution shows that, to leading order, the core scalar field, , exhibits a damped, wave-like variation in magnitude.

## 7. Summary

Three stochastic-based methods have been presented for solving unsteady scalar transport problems on bounded domains. The results, grouped according to Method number, are as follows:

A standard (Friedman 1975), backward-in-time solution to a Fokker–Planck equation is adapted to scalar transport problems, subject to Dirichlet conditions. The method, which has an intuitive geometric basis and been well studied within the mathematical community, has not received much attention by engineers and applied scientists. This paper attempts to make the method more accessible to these latter groups.

Three techniques are proposed for treating the practically important problem in which Dirichlet boundary conditions are

*a priori*unknown on a portion of the problem boundary. The most rigorous of these rests on calculation of the exit probability for RWs leaving the solution domain through the uncharacterized boundary.A simple method, based on quasilinearization, is outlined for tackling nonlinear transport problems, subject to Dirichlet conditions.

Grigoriu & Samorodnitsky (2003) method for solving steady mixed boundary value problems is extended to allow local, i.e. pointwise, solution of unsteady mixed problems. The method, which uses a simple reflection process for RWs impinging on Neumann boundaries, offers a straightforward alternative to existing local time-based approaches. Equally important, it is shown that Method II is equivalent to Milshtein's low-order scheme (Milshtein 1997) for mixed IBVPs; this equivalence provides the method with a rigorous theoretical basis.

Method II is used to derive a probabilistic-based analytical solution to Neumann problems; the solution incorporates the exit probability for RWs impinging on these boundaries. In two limits, using appropriate values for the reflection distance,

*ϵl*_{0}, the solution exactly matches continuum analytical solutions. Practically speaking, the technique offers a new and apparently novel approach for tackling both Neumann and Dirichlet problems (where application to the latter class of problems will be reported in a separate publication).The above comparison of probabilistic and continuum analytical solutions leads to consistent expressions for the off-set distance,

*ϵl*_{0}, used to reflect impinging RWs from Neumann boundaries. In the cases examined,*ϵl*_{0}is shown to be proportional to an incremental diffusion distance, , where 2*α*is the scalar diffusivity and Δ*θ*is the computational time-step size.A second, whole-field solution to mixed IBVPs is proposed. The method combines boundary Taylor expansions with Method I to arrive at a system of equations in the set of unknown scalar magnitudes and unknown surface fluxes extant, respectively, on Neumann and Dirichlet boundaries. The resulting linear system can then be solved via standard numerical approaches.

By choosing the Taylor expansion distance, Δ

*n*, to be of the order of*ϵl*_{0}, it is shown that Method III and a global version of Method II are computationally equivalent.

## Footnotes

- Received May 15, 2006.
- Accepted August 9, 2006.

- © 2006 The Royal Society