Interfacial mass transport in steady three-dimensional flows in microchannels

Rates of interfacial mass and heat transfer from uniaxial laminar flows through ducts at low Reynolds number follow the behaviour elucidated by Grætz: the local Sherwood (or Nusselt) number falls to an asymptotic value that is defined by the geometry and is independent of the Péclet number. In a previous study, we showed that analogous but distinct behaviour occurs in duct flows that exhibit Lagrangian chaos in the cross section due to spatially varying secondary flows: the Sherwood number follows the identical decay in the entrance region before passing to a Péclet-dependent value at larger axial distances; we called this behaviour, ‘modified Grætz'. Here, we investigate the generality of this behaviour in chaotic and non-chaotic flows for transfer to moving interfaces and across internal interfaces between convectively disconnected sets in the flow. We present theoretical predictions of the transfer rates for these cases and verify the accuracy of these predictions with a tracer-based, numerical simulation. We also present a theoretical criterion for the modified Grætz behaviour in terms of the axial length associated with mixing and the axial length associated with the return of depleted fluid to the reactive interface; we exploit simulation to verify this criterion.


Introduction
Interfacial transport is essential to the success of many processes that operate under laminar flow conditions. Heat and mass transfer to solid-liquid interfaces and across liquid-liquid interfaces are fundamental to heat exchanger design (Acharya et al 1992, 2001, Mokrani et al 1997, Peerhossaini et al 1993, electrochemical systems for analysis and energy production (Cohen et al 2005, Ferrigno et al 2002, Shrivastava et al 2008, separations with membranes (Shrivastava et al 2008) and without membranes (Aref andJones 1989, Brody andYager 1997), fabrication at fluid-fluid interfaces (Kenis et al 1999, and sensors involving interfacial reactions (Foley et al 2008, Golden et al 2007, Kamholz et al 1999, Squires et al 2008, Teles and Fonseca 2008, Vijayendran et al 2003. Many of these systems rely on the suppression of turbulence as a means of controlling fluid flow, whereas others are relegated to laminar flow by other constraints. These systems necessarily operate at low Reynolds number, Re = U H /ν < 100, but can nonetheless also operate at high Péclet number, Pe = UH /D > 100, where U is the average velocity of the fluid, H is the characteristic dimension of the system, ν is the viscosity of the fluid, and D is the diffusivity of the solute in question (or the thermal diffusivity of the fluid). While low Re leads to the laminar conditions that are either desired or otherwise required in these systems, high Pe generally leads to low efficiency (e.g. fractional conversion of reactant) in heat and mass transfer due to slow diffusive transport relative to the overall convective transport of the fluid. As Péclet numbers tend to be higher for mass transfer than for heat transfer in liquid flows, this work focuses on mass transfer processes; the results are nonetheless relevant to heat transfer as well.
The goal of this work is to understand how the character of a given laminar flow impacts rates of interfacial mass transfer, and to illuminate the design of flows that allow for high rates of interfacial transfer. Many laminar flows are uniaxial in character, although 3 three-dimensional flow can be induced through inertial effects at moderate Re (as in Dean's flows), by patterning of the bounding surfaces of the flow , or by driving the surfaces bounding the flow (as in a lid-driven cavity (LDC)). Parameters that influence the rate of transfer include: shear rates in the flow and particularly at the interface at which the transfer process is occurring; strength of the transverse components of a three-dimensional flow relative to the axial component; boundary conditions on the velocity at the interface in question; and the presence or absence of modulation or variation of the transverse flow along the axial dimension of the system. Metrics for evaluating the connection between the character of a flow and its interfacial transfer characteristics include: the mass transfer coefficient for transfer to the interface in question; the efficiency with which the flow homogenizes concentration differences in the bulk; and the structure and size of separate chaotic and non-chaotic sets that arise due to the flow.
We take as our system the class of steady flows in microfluidic channels of rectangular cross section (figure 1). These flows are sufficiently simple so as to be generalizable to the greater field of laminar flows, but also allow for significant complexity through the application of velocity boundary conditions along the perimeter of the channel. By applying nonzero conditions to the tangential components of velocity along regions of the bounding surfaces, recirculation in the cross section can be induced, giving the flow a three-dimensional character. This is the so-called LDC system, which can be applied in two or three dimensions. Three examples of three-dimensional LDC flows with mean axial velocity are shown in figure 1. The axial component of velocity is driven from left to right by a global difference in pressure, whereas the transverse components of the velocity are driven by slipping boundaries at the floor in figures 1(a) and (b), and at the floor and ceiling in figure 1(c). In figures 1(a) and (c), the transverse velocity field is modulated periodically along the length of the channel (with period 2L 1/2 , where L 1/2 is the half-cycle length). This modulation of an asymmetric flow in space or time gives rise to chaotic trajectories in the flow, enhancing mixing (Aref 1984). Such LDC systems and the chaotic properties of the flows that they produce have been studied extensively (Aidun et al 1991, Chien et al 1986, Leong and Ottino 1989. The chaotic nature of the flows can be tuned through variation of the parameters in figure 1: the degree of asymmetry, r, the half-cycle length, L 1/2 , and the transverse velocity at the moving boundary, u trans . The flow in figure 1(b) is three-dimensional, but without modulation of the flow in space or time; it is necessarily nonchaotic (Ottino 1989). We note that the red streamlines shown in the cross sections of the flows in figure 1 are the bounding lines of the vortices produced by these flows, connecting stagnation points of the transverse flow. While there is no convective transport across these streamlines, the modulation of the flows in figures 1(a) and (c) cause a periodic exchange of fluid between the left and right vortices. The red vertical streamline in figure 1(b) and the horizontal centreline in figure 1(c) are interfaces between vortices that are never modulated and therefore there is never convective transport across these streamlines; they are purely diffusive interfaces.
The LDC geometries used in this work generate bulk flows that are similar to the flows generated by the staggered herringbone mixer (SHM), a laminar, chaotic micromixer (Stroock et al 2002a). Patterning the surface of a microchannel with obliquely oriented grooves can impart a transverse component to the velocity near the grooved surface (Stroock et al 2002b). The net effect of these grooves on the flow in the bulk, to leading order, is equivalent to the effect of a slipping surface, as in the LDC. In this respect, we expect results for mass transfer to interfaces far from the driving surfaces in the LDC flows (such as the blue surface in figures 1(a) and (b) and the horizontal centreline in figure 1(c)) to reflect the effect of the transverse flows Characteristic parameters of the mixer are: channel height, H, and width, W, degree of asymmetry in moving boundary, r, half-cycle length, L 1/2 and transverse velocity at the driving boundary, u trans . (a) Chaotic transverse flow (r = 1/3) driven from the floor only; the reactive interfaces considered in this flow include the blue (stationary) interface and the green (moving) interface. In this flow, the degree of asymmetry is switched from r = 1/3 to r = 2/3 at the end of each odd half-cycle, and switched back to r = 1/3 at the end of each even half-cycle. (b) Non-chaotic transverse flow (r = 1/2) driven from the floor only; the reactive interfaces considered in this flow include the blue (stationary) interface and the green (moving) interface. In this flow, the transverse flow is not changed as a function of axial distance. (c) Chaotic transverse flow (r = 1/3) driven from the ceiling and floor in a mirror-symmetric manner; due to the symmetry of the driving velocities, the flow is symmetric about the centre plane at y = 0, and simultaneous transport to the bottom green (moving) interface and across the horizontal centreline is considered. generated by SHM-like devices, as illustrated previously (Stroock and McGraw 2004). We do not expect results for mass transfer to the moving boundaries in figure 1 to accurately reflect the effect on mass transfer to the grooved surface in SHM-like devices. The SHM grooves are composed of a collection of stationary surfaces whereas the moving boundary in the LDC is a single moving surface. While results for mass transfer to the moving interface cannot be directly related to those in the SHM, they are still instructive as a different type of interface for mass transfer and they may still shed light on some of the behaviours seen near the SHM grooves. Results for mass transfer across the diffusive interface at the horizontal centreline in figure 1(c), being far from the surfaces driving the transverse flow, should be indicative of the effects on mass transfer across any such interface in a symmetrically driven flow, whether it be driven by moving boundaries, SHM-like grooves, or any other convective process, such as an adaptation of the T-sensor system for separations (Brody and Yager 1997) or analysis (Kamholz et al 1999).
The problem of mass transport to the interfaces bounding the steady, pressure-driven flow in a rectangular duct bears a significant resemblance to the classic Graetz problem for heat transfer to the interfaces bounding the laminar, fully developed flow in a round pipe or between infinite flat pipes (Graetz 1885). By the Chilton-Colburn analogy between heat transfer and mass transfer at high Pe, we can map the Graetz result for heat transfer between a fluid with uniform inlet temperature and the solid bounding interface at some other uniform surface temperature onto the problem of mass transfer between a fluid with uniform inlet solute concentration and the solid bounding interface at some other uniform surface concentration (Bird et al 1960). We will therefore consider the classic Graetz problem for the flow between infinite flat plates using the language of mass transfer to consider an instantaneous chemical reaction at the stationary interface resulting in zero concentration at the interface.
As the fluid at uniform inlet concentration encounters the reactive interface, there is a flux of solute to the interface from the nearby fluid due to the difference in concentration between the fluid and the interface. The fluid adjacent to the interface becomes increasingly depleted of solute, and a concentration boundary layer forms and grows in the direction perpendicular to the interface with increasing axial distance. The flux to the reactive interface decreases as the growing boundary layer causes a decrease in the gradient of concentration at the interface. A mass transfer coefficient relating the flux at the interface to the concentration driving force is defined by: where k(z) is the mass transfer coefficient, J(z) is the flux of solute to the interface, C surface is the surface concentration at the reactive interface (taken to be zero), and C bulk (z) is the cup-mixing average concentration: The mass transfer coefficient can be non-dimensionalized for comparison across systems of different size scales and solute diffusivity, and can be written by conservation of mass of solute: where Sh(z) is called the Sherwood number for mass transfer, H is the separation between the plates and Z = z/PeH is the scaled axial distance. For short axial distances, the flux to the 6 reactive interface decreases due to a decrease in the Sherwood number, while the concentration driving force stays nearly constant. This region is called the entrance region. As the thickness of the concentration boundary layer approaches the separation between the plates, however, the concentration profile develops an asymptotic shape and the Sherwood number becomes constant while the average concentration decreases with axial distance (see figure 2). This region is called the asymptotic region and is characterized by an asymptotic value of the Sherwood number, Sh ∞ , and the axial distance at which the transition occurs is called the entrance length, L ent . For uniaxial flow between infinite flat plates, L ent /PeH and Sh ∞ are universal constants. For uniaxial flow in a rectangular duct, L ent /PeH and Sh ∞ have additional geometry dependence, but are independent of Pe for a given geometry.
In the presence of a three-dimensional flow, the mass transfer problem considered above takes on new characteristics that can lead to larger values of the Sherwood number and higher flux of solute to the reactive interface. With flows like those created in the geometries shown in figure 1, the depleted fluid near the reactive interface can be swept away from the interface and into the bulk. If the bulk is also sufficiently well mixed, the depleted fluid will be homogenized with the bulk before it returns to the surface. In this way, solute is consistently brought near to the reactive interface by convection rather than by diffusion across a thick boundary layer. The transverse flow stops the growth of the boundary layer due to the movement of depleted fluid away from the reactive interface: once the fluid has been swept across the width of the reactive interface, it detaches and is incorporated into the bulk stream (see figure 2). The entrance region in such a flow is similar to that in a uniaxial flow, but the asymptotic state is reached after a shorter scaled entrance length, L ent /PeH , and at a higher asymptotic value of the Sherwood number, Sh ∞ ; L ent /PeH and Sh ∞ are both functions of Pe. We call this effect the modified Graetz behaviour because it shares the hallmarks of the classic Graetz solution, but it introduces dependence on Pe into L ent /PeH and Sh ∞ ; these characteristics allow for significant increases in total reactive flux to the interface and efficiency of usage of the solution of reactants.
In our previous work, we demonstrated the existence of the modified Graetz behaviour in the context of transport to the stationary interface far from the surface that drives the flow, as in figure 2 (Kirtland et al 2006). The present work has two main purposes: to extend the investigation of the modified Graetz behaviour to include internal diffusive interfaces and moving reactive interfaces; and to develop a more complete understanding of the origin and generality of the modified Graetz behaviour. The first goal (extension to new interfaces) includes the development of theoretical predictions for the effect of transverse flows on mass transfer to or across the interfaces in question, and then the validation of those predictions from numerical simulations. The second goal (understanding the origins and generality of the modified Graetz behaviour) involves the assessment of the dominant length scales in a given flow, and the determination of whether the processes that drive mixing in the bulk are sufficiently fast to keep the processes that drive depletion of the boundary layer from reinforcing and decreasing the efficiency of the system.
In order to place this work in the greater context of heat and mass transfer in laminar flow, we now consider the existing work in the field. Theoretical developments are discussed first, and experimental and numerical results follow. As discussed above, the classic problem of heat transfer from a fluid undergoing fully developed laminar flow to the surface of a round pipe was solved by Graetz in 1885. Further insight was added by Lévêque in 1928 in the form of a similarity solution for the temperature field in a shear flow, which is a good approximation of the flow near the pipe surface, and therefore is meaningful for thin boundary layers, i.e. in the entrance region. These ideas were mapped onto the problem of transfer across a fluid-fluid interface at the centreline of the flow in a rectangular duct by Ismagilov et al (2000), demonstrating that the rate of mass transfer across the interface, and therefore the rate of homogeneous reaction of initially segregated reactants, scales with the square root of Pe in the plug-like flow in the centre of the channel, but scales with the cube root of Pe near 8 the stationary duct walls. Chang and Sen (1994) called attention to the types of resistances to transfer from flows to solid surfaces and across fluid-fluid interfaces, identifying boundary layers, recirculation regions and stagnation streamlines, and also noted that chaotic mixing leads to lower resistances than non-chaotic mixing. As the introduction of chaotic trajectories requires modulation, and, therefore, a frequency of modulation, many groups have investigated the optimal frequency for increasing rates of transfer in a given system (Bryden and Brenner 1996, Ghosh et al 1992, Horner et al 2002. Although many have noted significant increases in rates of transfer with the introduction of chaotic trajectories (Bryden and Brenner 1999, Ganesan et al 1997, Lefevre et al 2003, others (ourselves included) have cautioned against equating increases in mixing efficiency due to chaotic flows with increased rates of transfer (Ganesan et al 1997, Ghosh et al 1998, Kirtland et al 2006. Experimentally, many investigators have harnessed the controllability of laminar flow to their advantage in systems depending on interfacial heat and mass transfer: liquid-liquid extraction (Brody and Yager 1997); homogeneous reaction of initially separate reactant streams , Kamholz et al 1999; microfabrication at fluid-fluid interfaces (Kenis et al 1999 and laminar flowbased fuel cells (Cohen et al 2005, Ferrigno et al 2002. All of these systems, however, suffer some loss of efficiency due to depleted regions near the interfaces at which transfer occurs.
Experimental efforts to understand the influence of chaos on interfacial transfer are many and varied in their conclusions. Researchers working with Dean's flows in coiled tubes have cited enhanced heat transfer due to chaos (Peerhossaini et al 1993), although the increases in rates of transfer are relatively small in magnitude (13-27%, Mokrani et al 1997) and in scaling with the Prandtl number, Pr = ν/α, where ν is kinematic viscosity and α is thermal diffusivity (Pr 0.09 , Acharya et al 2001). For transport to the surface of confocal ellipses, a 100% increase over conduction is observed by Saatdjian et al (1996). For mass transfer to the surfaces of a sinusoidal wavy-walled channel, chaos leads to higher rates of transfer, and the effects are associated with the flow separation and oscillation induced by the wall geometry (Nishimura 1995, Nishimura and Kojima 1995, Nishimura et al 1993. Finally, Shrivastava et al (2008) were able to predict to within 4-12% the enhancements seen in experiments with ultrafiltration membrane spacers using a simple model, even though this model assumes instantaneous mixing of the fluid. While it is difficult to incorporate the many experimental, numerical and theoretical results available in the literature into a single conclusion, a theme emerges: chaos can lead to significant increases in transfer to moving interfaces and internal, fluid-fluid interfaces (Bryden and Brenner 1999, Ganesan et al 1997, Lefevre et al 2003, while for transfer to solid, stationary interfaces, chaos leads to little or no increase in rates relative to analogous, nonchaotic flows (Acharya et al 2001, Ghosh et al 1998, Kirtland et al 2006, Mokrani et al 1997, Peerhossaini et al 1993.
In order to develop a better understanding for the mechanism behind the increases in rates of transfer that are seen in the presence of transverse flow, we use simulation code designed to track particles in a three-dimensional flow field to simulate the convection-diffusion-reaction process at stationary and moving solid interfaces, and also to simulate the convection-diffusion process at an internal, fluid-fluid interface. We look for signs of the modified Graetz behaviour described above in the results for mass transfer to these interfaces, and, perhaps more importantly, we look for situations in which the modified Graetz behaviour breaks down at large axial distance and at high Pe. In relating the success or failure of the modified Graetz behaviour to the character of the flow, we propose a more complete picture of the mechanism behind the modified Graetz behaviour and behind the maintenance of high rates of mass transfer, in general.

9
The paper is structured as follows: in section 2, we present a definition of the standard Graetz and modified Graetz behaviours described above and review the results from our previous work (Kirtland et al 2006). We then proceed to develop theoretical results for the extension of the modified Graetz behaviour to moving reactive interfaces and internal diffusive interfaces, as well as a theory for the conditions on the existence of the modified Graetz behaviour in terms of characteristic length scales in the flow. Section 3 describes the numerical methods used to test theoretical predictions. In section 4, we present and discuss results for interfacial mass transfer in the flows in figure 1, as well as results for characteristic length scales in the flow, and an extension of the modified Graetz framework to a broader class of flows. Finally, we present conclusions and discuss the applicability of our analysis to the field of interfacial transfer phenomena more generally.

Theory
In this section, we present the theoretical framework in which we will interpret our results in section 4. Starting with the governing equation for the convection-diffusion process, we derive forms for the local value of the Sherwood number for transport to the blue and green interfaces in figure 1 and for transport across the plane of symmetry in figure 1(c). We also develop correlations for the asymptotic value of the Sherwood number that are valid beyond the entrance length, once the concentration profile has converged to its asymptotic shape. Finally, we develop an argument for the criteria for the existence of modified Graetz behaviour at arbitrarily large axial distance and arbitrarily large Pe.

Mass transfer to stationary interfaces in LDC flows
As discussed in section 1, the groundwork for the understanding of transport to stationary solid interfaces in fully developed uniaxial flow was laid by Graetz (1885). The convection-diffusion equation is the governing equation for the process (axes defined in figure 1): where u is the velocity, C is the concentration, and u z is the axial component of the velocity. Symmetry is assumed in the x-direction and streamwise diffusion is assumed to be negligible compared with spanwise diffusion, which is a reasonable assumption at high Pe. A solution of the governing equation by separation of variables is possible for the case of uniform inlet concentration and uniform concentration at the interface, resulting in an exponential function of z and a power series solution in y. The full solution includes a sum over eigenmodes, but at large axial distance, only the mode with the smallest eigenvalue persists. This mode determines the shape of the asymptotic state mentioned in section 1, and the asymptotic value of the Sherwood number. For a thin boundary layer near the interface, as is the case in the entrance region, transfer occurs in a shear flow with zero velocity at the interface. In 1928, Lévêque solved the problem of transfer from a fluid undergoing shear flow to the stationary surface bounding the flow by use of a similarity variable, making it simple to calculate the Sherwood number in the entrance region as a function of axial distance: where Sh ∞ and L ent /PeH are geometry dependent but independent of Pe, with Sh ∞ ≈ 2.55, L ent /PeH ≈ 0.06 and B 0 ≈ 1 for the uniaxial flows considered in this work.
In systems with transverse as well as axial components of velocity, there exists the possibility of convecting the depleted fluid near the reactive interface into the bulk, such that the boundary layer does not become successively thicker with increasing axial distance. In order for this modified Graetz behaviour to hold, the transverse flow must be strong enough relative to the axial flow to sweep fluid away from the reactive interface over an axial distance that is short compared to the entrance length for the standard Graetz behaviour, otherwise, the effect of the transverse flow will be obscured by the rapid growth of the standard boundary layer. In order for this behaviour to continue to hold for large axial distances, the transverse flow must also homogenize the depleted fluid with the bulk fluid before it returns to the reactive interface. If these conditions are satisfied, the asymptotic state of the concentration distribution takes on a different form determined by the transverse flow. The Sherwood number in the entrance region is identical to that in the uniaxial case, but the entrance length is shorter and independent of Pe: L ent ∼ Wγ axial /γ trans for width of interface, W , and axial and transverse shear rates at the interface,γ axial andγ trans . The asymptotic state involves a concentration profile that is selfsimilar (periodically self-similar for modulated flows) with a mean concentration that continues to decay as in the uniaxial asymptotic region, but the asymptotic value of the Sherwood number for transfer to the stationary surface, Sh ∞ stat , becomes Pe-dependent: where Pe trans = Peu trans /U . This result is the essence of the modified Graetz behaviour as presented in our previous work (Kirtland et al 2006). We proceed now to extend the modified Graetz result to other interfaces in the flows considered in this work.

Mass transfer to non-stationary interfaces in LDC flows
2.2.1. Transport to diffusive internal interfaces. In order to extend the above development to moving interfaces, we start with the simple example of diffusion across the horizontal interface between two streams of different initial concentration in a uniaxial plug flow in the z-direction. We can approach this interface in much the same way that we approached the simple shear flow as a model of the flow near a generalized solid stationary interface (Kirtland et al 2006). For transfer across the x-z plane (analogous to the horizontal centreline in figure 1(c) in a uniaxial plug flow), the problem is symmetric in x. The transverse components of the velocity are zero, and we assume that at large Pe we can neglect streamwise diffusion. The convection-diffusion equation is solved by similarity solution to give: where u z is the axial velocity, D the diffusivity and C 0 the initial difference in concentration across the centreline. The form for the local Sherwood number Sh(z) is where Pe plug = u z H/D. The local mass transfer coefficient decreases with increasing axial distance due to a thickening boundary layer at the interface, which results in a decreased gradient of concentration driving the transfer process. By the same logic, we can consider transport to the diffusive internal interface (y = 0) in the flow in figure 1(c). It should be noted, however, that (8) is the result for a true plug flow, and the extension of this result to the diffusive interface in figure 1(c) is only approximate due to the mixed scaling found in the vicinity of the sidewalls, as demonstrated by Ismagilov et al (2000). There is a plane of symmetry at the horizontal centreline of the channel due to the mirror symmetry of the driving velocities at the ceiling and floor of the channel. The velocity normal to this plane must be zero, and the derivative of the velocity tangent to the plane must be zero. Zero velocity normal to the plane results in an interface across which solute only diffuses; there is no convective transport across this centre plane. The interface acts like the membrane in a coflowing membrane separator (Brody and Yager 1997), and such interfaces act like a 'virtual membrane' in the context of fluid-based fuel cells (Cohen et al 2005, Ferrigno et al 2002.
The vanishing derivative of the tangential velocity gives us effectively a plug flow in both the axial and transverse directions near the centreline. We can consider the growth of a boundary layer in the axial and transverse directions. The governing equation is again the convection-diffusion equation, but now there are two convective terms: where u centre is the transverse velocity at the centreline due to the driving velocity at the floor and ceiling, u trans , and u axial is the axial velocity at the centreline, which is roughly twice the average axial velocity U. Non-dimensionalizing (9) where C o is the initial concentration difference across the interface, W is the channel width, H is the channel height, L is the characteristic axial length over which the boundary layer develops due to the axial flow (L = L ent ) and Pe = UH /D. Transverse convection will dominate development of the boundary layer if the axial convective term is small compared with the transverse convective term. This situation occurs if L ent /W 2U/u centre . From past results (Kirtland et al 2006), we expect the axial distance required to develop a boundary layer through a height h in the channel by diffusion in a uniaxial flow to be roughly L ent = 0.06Peh, and in this case, h is the half-height of the channel: h = H/2. Therefore we expect the transverse convection term to dominate for,

12
where Pe centre = Peu centre /U . In this limit, we can neglect the axial convection term in (9), such that the governing equation becomes As for (7) above, there is a similarity solution for (12) The form of Sh(x) developing across the channel is now Physically, the concentration profile and the corresponding Sh centre (x) continue to develop across the channel until the fluid is swept off the internal interface and back into the bulk. This departure from the interface arrests the growth of the depleted region and leads to the plateau value of the Sherwood number, in a manner consistent with the modified Graetz behaviour (figure 2). This plateau value is equal to the average of Sh centre (x) over the width of the channel.
As can be seen in the streamlines in figure 1(c), the effect of the smaller vortex at the centreline is minimal; the larger vortex accounts for most of the centreline velocity. Averaging Sh centre (x) over the width of the channel, we have the following correlation for the asymptotic value of the Sherwood number at this symmetry plane: This asymptotic value of the Sherwood number holds once the boundary layer near the internal interface has assumed its asymptotic shape, after the end of the entrance region, with length L ent ∼W 2U/u centre . Within the entrance region, Sh centre (z) takes on the same values as the local Sherwood number in a plug flow as in (8) with Pe plug = u z H/D = 2UH /D = 2Pe.

Transport to moving interfaces.
We now consider transport to the moving interfaces in the flows in figures 1(a)-(c). The moving boundary that drives the transverse flow is unusual in that it is a mixed slip, stationary surface. In the axial direction it is stationary and presents roughly a simple shear flow near the surface, but in the transverse direction, it moves with uniform velocity. This motion creates a transverse shear flow in its vicinity: a shear flow with nonzero velocity at the wall and a negative shear rate with respect to distance from the wall. This flow can be approximated as u x (y) ≈ u trans −γ trans y near the wall, where y is the distance from the wall. This shear flow term significantly complicates the solution of the governing equation, as we now have a sum of plug-like and shear-like components. If we consider the flow very near the moving boundary, say within a distance y u trans /γ trans , we can assume that the plug-like component dominates and we can drop the shear-like term. We will consider the validity of this assumption by finding the approximate maximum thickness of the depleted region as compared with this condition on y .
Proceeding with the plug flow solution from above, we take the form for Sh(x) from (14), replacing Pe centre with Pe trans = Peu trans /U and noting that this interface has only half the resistance to transfer compared with the centreline case above: in the centreline case, solute must cross twice the boundary layer thickness from the bulk on one side of the interface to the bulk on the other side of the interface, and therefore the form of Sh(x) is twice as large at the moving reactive interface as it is at the diffusive internal interface: In this case, we average across the half-width, because the smaller vortex is strong enough to contribute to the development of the depletion layer ( figure 1(a)). We find Again we arrive at a prediction for the prefactor of Sh ∞ trans and for its scaling with Pe and u trans . In order to consider the validity of our assumption above about being sufficiently near the wall, we consider the approximate thickness of the depleted region, δ(x). For a reasonable estimate of this thickness, we take: In the flows considered in this study and depicted in figure 1, the shear rate at the moving wall is approximately 5u trans /H . To ensure that the boundary layer does not grow thick enough to sample a region where the shear term is greater than one-tenth the size of the plug flow term, we need the following: Pe trans > 50 2 π 2 W H ∼ 4000.
For Pe trans larger than about 4000, the region of the flow sampled by the boundary layer sees a nearly uniform plug flow with a deviation of less than 10%, and we will expect our simple plug flow arguments to hold for Pe trans > 4000 and that (17) should provide a reasonable estimate of Sh ∞ trans .

Summary of modified Graetz correlations.
We summarize the correlations presented here and elsewhere (Kirtland et al 2006) for the asymptotic Sherwood number due to the modified Graetz behaviour. The asymptotic values in table 1 are expected to hold at axial distances beyond the entrance length. Within the entrance region, all interfaces that are stationary with respect to z have the same form for Sh(z): and interfaces that are plug-like in z have a similar form for Sh(z): This allows for an estimate of the approximate entrance length by setting the Sherwood number at the transition from the entrance region equal to the value in the asymptotic region:

Requirements for the modified Graetz behaviour
With the predictions derived in the previous section for flows that have the potential to display the modified Graetz behaviour, we now consider the following questions: for a flow that demonstrates the hallmarks of the modified Graetz behaviour (departure of Sh(Z) from uniaxial result at small to intermediate axial distance) what is required for the modified Graetz behaviour to hold for arbitrary axial length? For a flow that demonstrates modified Graetz behaviour at some Pe, what is required for the modified Graetz behaviour to be valid for arbitrarily large Pe?
To address these questions, we compare characteristic length scales in the flow. We start with the assumption that the flow in question has the ability to sweep the depleted fluid near the reactive interface into the bulk. If this fluid is immediately brought back to the leading edge of the reactive interface, the depletion will be reinforced and the thickness of the boundary layer will continue to grow, as in the uniaxial case. If, however, the depleted fluid remains in the bulk long enough before returning to the reactive interface that it homogenizes with the bulk, there will be no memory of the original depletion; the fluid impinging on the reactive interface will be at the local, cup-averaged concentration of the bulk. We therefore define two characteristic axial length scales: the mixing length, L mix (Pe), and the return length, L return (Pe). At a given Pe, if the flow displays modified Graetz behaviour at small to intermediate axial lengths and L mix is shorter than L return at that Pe, the depleted fluid will always be reconcentrated before it returns to the reactive interface and the modified Graetz behaviour will hold for arbitrarily large axial distances. If the modified Graetz behaviour holds for arbitrary axial length at some Pe and L mix scales more slowly than L return with Pe, then the modified Graetz behaviour will hold for arbitrarily large Pe. We now define the mixing length and return length more fully and consider how they are expected to depend on Pe.
2.3.1. Mixing length: definition and scaling. For a diffusive scalar with a non-uniform inlet distribution, diffusion will tend to decrease the heterogeneity of concentration through flux of solute along the concentration gradients in the cross section, making the fluid more homogeneous or mixed. As a measure of the degree of mixing, we take the velocity-weighted variance of concentration in the cross section: where the cross section has been discretized into N bins, C i (z) is the concentration in bin i, u i is the average velocity in bin i, A i is the area of bin i, U is the average velocity for the whole cross section, C cup is the velocity-weighted average concentration for the whole cross section (as in (2)), C max is the maximum inlet concentration, and A is the cross-sectional area. This measure goes from an initial value of 1 at the inlet to 0 at infinite length (complete mixing).
We employ a velocity-weighted measure of mixing as opposed to the usual uniformly weighted variance because, if the fluid at a particular axial distance were sampled for some period of time, the deviations from the average concentration would contribute to the overall variance of concentration in proportion to the local axial velocity. For much the same reason, we chose the 'cup-mixing' average concentration in (2). In order to identify the scaling of the mixing length with Pe, we plot the variance from numerical simulation versus scaled axial distance, z/ f (Pe), for various forms of the scaling function, f (Pe); we take the f (Pe) that provides the best, qualitative collapse of σ 2 (z/ f (Pe)) for multiple values of Pe to represent the Pe-scaling of L mix . For a uniaxial flow, mixing occurs purely by diffusion and the mixing length is expected to scale linearly with Pe. For an efficiently mixed flow such as the chaotic flow in figure 1(a), the mixing length is expected to scale as either ln(Pe) (Ranz 1979, Villermaux and Duplat 2003, Villermaux et al 2008 or as a weak power law in Pe (Simonnet and Groisman 2005). The non-chaotic three-dimensional flow in figure 1(b) is expected to show scaling with Pe that falls somewhere between these two extremes.

Return length: definition and scaling.
In determining the dependence of the return length L return on Pe, we consider the trajectory of a particle of solute that has reacted at the reactive interface. We call this particle a depleton, as it represents a hole in the concentration distribution left by the depletion of fluid near the interface. The trajectory of this depleton will eventually return to the reactive interface, at which point the absence of an actual reactive particle leads to a decrease in the gradient of concentration at the interface. The return length for a depleton is therefore defined as the axial distance between consecutive collisions with the reactive interface (see figure 3). Due to the diffusive displacements of particles of solute in the flow and the non-uniformity of the three-dimensional flows considered, there is not a single, well-defined return length, but rather a distribution, as illustrated by the sample trajectories in figure 3. There are some depleton trajectories that have a very short return length: these trajectories are associated with the initial depletion of the fluid as it passes near the reactive interface, i.e. the growth of the boundary layer, and they never actually enter the bulk. Longer return lengths correspond to depletons that have been swept into the bulk by the transverse flow, and subsequently returned to the region near the reactive interface by the transverse flow. These longer return lengths represent the part of the distribution that is relevant to the question at hand, as they correspond to the trajectories that have a chance of being mixed with the bulk before they return.
We expect L return to show the same scaling with Pe as the axial length associated with decreasing the bulk concentration, because the transfer process that takes the depleton into the bulk and returns it to the reactive interface is the same convective process that brings unreacted solute particles from the bulk to the reactive interface. The expected scaling for a decrease in the average concentration by a factor of e is derived from the dependence of the concentration on axial distance and on the average Sherwood number for transport to the interface: where L 1/e is the length required to decrease the concentration by a factor of e and Sh 1/e is the Sherwood number averaged axially from the inlet up to length L 1/e , which, in the case of a short entrance region, can be taken as the asymptotic value of Sh(z), Sh ∞ , for transport to the interface. If Sh ∞ ∼ Pe 1/3 , as in the case of transport to a stationary interface, the return length should scale like L return ∼ Pe 2/3 . If Sh ∞ ∼ Pe 1/2 , as in the case of transport to a moving or diffusive internal interface, the return length should scale like L return ∼ Pe 1/2 .

Model and simulation details
Simulations of mass transport to interfaces were carried out by advecting diffusive tracer particles in an analytical flow field that is the approximate solution for the pressure-driven Stokes flow in a three-dimensional LDC. The basis of these simulations is identical to that which was presented in our previous work (Kirtland et al 2006). The velocity field is decoupled into axial and transverse components, where the axial component is the familiar Poiseuille solution for axial flow in a duct: where γ n = (2n − 1)π/2, and the transverse components are calculated by solving the biharmonic equation for the streamfunction in the cross section: The specifics of the velocity calculations are laid out in detail elsewhere (Stroock and McGraw 2004).
For the measurements of concentration and interfacial flux, particles were seeded into the flow at the inlet and advected through a distance equivalent to 100 cycles of the modulated lid driven cavity flow shown in figure 1. The initial positions were chosen to simulate the flux profile impinging on the inlet for a uniform inlet concentration in fully developed Poiseuille flow. For any given time step, the simulation integrated the velocity by use of a fifth-order Runge Kutta scheme with an adaptive time step. The velocities used in the Runge Kutta step were interpolated from a fine (512 × 512 or finer) grid of velocities calculated beforehand from the double Fourier expansion due to Stroock and McGraw (2004). A diffusive step was also employed by adding pseudorandom displacements in x, y and z chosen from Gaussian distributions with width (2D t) 1/2 , where D is the diffusivity of the particle and t is the current time step.
The locations of collisions with the walls were recorded, as well as the positions in the cross section as particles reached specific axial distances (usually the end of each cycle of the mixer or an equivalent distance in uniaxial or unmodulated flows). To simulate instantaneous kinetics, each collision was taken to represent the irreversible destruction of that tracer. This information allowed for the calculation of average concentration as a function of axial distance, C cup (z), which was proportional to the number of unreacted particles remaining at a given axial distance, and from that data we calculated the local Sherwood number as a function of axial distance as in (3).
We also calculated concentration profiles in the cross sections as a function of axial distance by binning the unreacted tracers in the cross section to find axial convective flux; concentration within each bin was then found by dividing this flux by the average axial velocity in the bin, calculated from (25). The cross-sectional positions were also used for the mixing length calculations. A modification was made to the code to record multiple collisions at a given wall: a distribution of return lengths was generated by recording the first 100 000 collisions of a tracer with each wall and analysing the distribution of distances between collisions; after each collision, we considered the particle to be a 'depleton' generated by that reaction event.

Transport to stationary, reactive interfaces
In this section, we present profiles of the concentration in the cross section of the channel and plots of the local Sherwood number, Sh(Z). We also present the concentration profiles as animations to demonstrate more effectively the evolution of the boundary layer structure. In figure 4 and movie 1 (available from stacks.iop.org/NJP/11/075028/mmedia), we show the evolution of the boundary layer for reaction at the blue (stationary) interface in the chaotically stirred flow shown in figure 1(a), with ten frames during each of the first ten cycles. The movie highlights the oscillating motion of the flow and how it sweeps the depleted fluid away from the reactive interface. The structure of the concentration boundary layer is well developed by the end of the first cycle, once the boundary layer has grown from the top left and right corners to the middle of the top interface. The transverse flow then sweeps the depleted fluid down into an oscillating tail and mixes it into the bulk.
The frames of all of the following animations (movies 2-4) correspond to the end of consecutive cycles in the modulated flows and corresponding axial distances in the unmodulated flows. In these movies, the concentration profiles in the modulated flows appear to be constant, despite the existence of variations throughout each cycle, as seen in movie 1. This effect is due to the fact that the profiles evolve in a periodic manner with the same period as the variations in the flow: C(x, y, z + 2n L 1/2 ) = C(x, y, z) f (n) where C(x, y, z) captures the periodic profile and f (n) captures the asymptotic decrease in the average concentration.
As illustrated in figure 5 and movie 2, more of the solute has reacted at the stationary interface at a given axial distance in the three-dimensional flows (second and third rows) than in the uniaxial flow (first row). The difference between the average concentration in the chaotic and non-chaotic flows is small relative to the difference due to the presence of a transverse flow. Figure 5 also shows that the transverse flows lead to an asymptotic shape of the boundary layer after a short axial distance (approximately 1 cycle = 1 frame in the movie), while the boundary layer in the uniaxial flow continues to grow throughout the length of the channel (100 frames). This is particularly clear in the third column of figure 5, where we have normalized the width-averaged concentration by the bulk average concentration: with the exception of the left-most pixel (representing the initial condition) the shape of the width-averaged concentration profile is constant through the entire length of the device for the cases with transverse flow. The short entrance lengths (L ent 2L 1/2 ) in the three-dimensional flows is also consistent with our previous observations (Kirtland et al 2006). Finally, we note that while the depleted region in the non-chaotic flow reconnects with the reactive boundary after circling the vortices, there is no sign of reconnection of the depleted region in the chaotic flow. The depleted fluid that is swept down through the bulk is twisted into the large vortex as the flow is modulated at the end of each half-cycle (see movie 1 at Pe = 10 4 to observe this process more clearly). This depleted region disappears into the bulk rather than returning directly to the reactive interface. Figure 6 presents the local Sherwood number, Sh(Z), as a function of scaled axial distance calculated from data such as that in figure 5 and movie 2. In the case of uniaxial flow (black symbols in figure 6), Sh(Z) falls on a single curve for all Pe, from 10 2 to 10 6 . This curve for Sh(Z) at all Pe has an entrance region that shows the expected scaling of Z −1/3 , an entrance length of L ent /PeH ∼0.06, and an asymptotic value of Sh ∞ ∼2.55, as in Kirtland et al (2006). In the case of a chaotic transverse flow (coloured symbols in figure 6(a)), Sh(Z) has a shape similar to that of the uniaxial case, with an entrance region where Sh(Z) decays and an asymptotic region where Sh(Z) becomes constant. For Pe < 10 5 , the entrance region is visible and falls on the same curve as the uniaxial case, and we expect that the entrance region of the cases with Pe > 10 5 lie on this curve as well. We also note distinct behaviour in the presence of transverse flow: firstly, the concentration distribution reaches its asymptotic state at shorter scaled axial distances, L ent /PeH as Pe increases; this short entrance length leads to a nearly constant value of Sh(Z > L ent /PeH ) = Sh ∞ stat through most of the length of the channel. Secondly, the value of Sh ∞ stat is Pe-dependent and can be much larger than in the uniaxial flow. The transverse flow maintains this higher rate of transfer by sweeping the depleted solution off of the reactive interface and thus maintaining a thinner depleted boundary layer, with a thickness that decreases with increasing Pe. The evenly spaced plateaus on this logarithmic plot imply a power law scaling of Sh ∞ stat with Pe, as will be confirmed in figure 11. Figure 6(b) presents Sh(Z) for the non-chaotic flow in the third row of figure 5 and movie 2. We see much the same effect due to transverse flow as was seen in figure 6(a): an entrance region 20 Figure 6. Local Sherwood number, Sh(Z), as a function of scaled axial distance, Z = z/PeH , for transport to the blue (stationary) interface in the flow shown in figures 1(a) and (b). Colours signify Pe (10 2 (purple), 10 3 (red), 10 4 (green), 10 5 (blue), 10 6 (grey)) at constant u trans = U/5, and the black dots signify the collection of all Pe for u trans = 0 (uniaxial Poiseuille flow). (a) Flow as in figure 1(a) with r = 1/3, u trans = U/5 and L 1/2 = 10H . (b) Flow as in figure 1(b) with r = 1/2, and u trans = U/5. that collapses to the uniaxial case (visible for Pe < 10 5 ) out to a scaled distance that is shorter than in the uniaxial case, followed by a Pe-dependent value of Sh(Z) beyond the entrance length L ent /PeH . For Pe < 10 5 , Sh(Z) reaches a well-defined plateau over two orders of magnitude in axial distance, with similar values to those seen for the chaotic flow ( figure 6(a)). At higher Pe, Sh(Z) does not plateau; rather, it continues to decrease with increasing axial distance. This decrease in Sh(Z) with increasing axial distance is due to the reinforcement of the boundary layer as depleted fluid returns to the reactive interface without being mixed with the bulk, as seen in the later frames of movie 2. For Pe < 10 5 , the depleted fluid is sufficiently reconcentrated as it circles the vortex that returns it to the surface, but for Pe > 10 5 , mixing from the cores of the vortices to the depleted fluid surrounding the cores is too slow to reconcentrate the fluid sufficiently. While a specific value of Sh ∞ stat cannot be defined in these cases due to the continuing decrease of Sh(Z), it is instructive to note that Sh(Z) is similar in magnitude for the chaotic and non-chaotic flows. As noted by others (Acharya et al 2001, Ghosh et al 1998, Kirtland et al 2006, Mokrani et al 1997, Peerhossaini et al 1993, from a globally averaged Sherwood number or from a total integrated flux over the entire reactive interface, it is very difficult to distinguish the transport process in figure 6(a) from that in figure 6(b). The non-chaotic flow fails to exhibit modified Graetz behaviour for Pe > 10 5 , but this failure is subtle in the sense that in many cases it would not even be measurable from an external, global measurement.

Transport to moving, reactive interfaces
In figure 7, the evolution of the concentration profile is depicted as in figure 5, but in this case solute reacts at the moving interface at the floor of the channel. The uniaxial flow presents the same behaviour that we saw in figure 5 because the uniaxial case is vertically symmetric, such that transfer to the top, stationary interface (as in figure 5) is equivalent to transfer to the bottom, stationary interface (as in figure 7). The chaotic flow keeps the boundary layer extremely thin (on the order of one pixel) and the rate of transport to the moving interface is so high that essentially all solute has reacted within the first half of the device. This high rate of transfer is due to both the magnitude and the plug-like character of the velocity at the interface. The non-chaotic flow looks quite similar initially, but as the depleted fluid encircles the vortices and returns to the reactive interface, the depletion is reinforced and the boundary layer continues to grow thicker with axial distance. A region of high concentration persists at the core of each vortex; these regions are not convectively connected to the depleted fluid outside of the core and thus cannot be efficiently mixed with the depleted fluid. Both the chaotic and non-chaotic flows send the newly depleted fluid to the region near the far wall, producing what is effectively an inverted boundary layer: the transverse flow moves the depleted solution away from the reactive interface to form a depleted zone near the far, non-reactive boundary; this zone is similar in shape to the usual boundary layer formed when solute reacts at the stationary interface. This inversion has the effect of keeping the concentrated fluid near the bottom reactive interface, thereby increasing transfer to the interface. The chaotic flow distinguishes itself by twisting the oscillating tail of depleted solution seen in figure 4 (of similar shape in figures 5 and 7) into the bulk rather than allowing it to reconnect to the surface as it does in the non-chaotic flow. The concentration profile in the chaotic flow again reaches an asymptotic shape within the first several cycles, as is evident in the normalized side view. In the non-chaotic flow, the boundary layer near the reactive interface and the inverted boundary layer near the ceiling both continue to grow thicker through the length of the device. The asymptotic state in the non-chaotic flow is reached only once the boundary layer near the reactive interface has grown thick enough to reach the core of the vortices, similar to the standard Graetz behaviour in the uniaxial case. Figure 8 presents Sh(Z) as a function of scaled axial distance for transport to the moving interface as shown in figure 7 and movie 3. The behaviour observed in figure 8(a) in the presence of chaotic transverse flow shares several important features with the results of figure 6(a) for transport to the stationary interface. The entrance region again lies on the curve for uniaxial flow, although, due to significantly higher rates of transfer at the moving interface, this collapse nonchaotically stirred with r = 1/2 and u trans = U/5. The first column presents the concentration, C(x, y, z), in the cross section at a distance z downstream that is equivalent to fifty cycles of the mixer (z = 100L 1/2 = 1000H ). The second column shows the width-averaged concentration in the cross section,C(y, z), as a function of axial distance, and the black line marks the position along the length of the channel corresponding to the cross-sectional profile in the first column. The third column is equivalent to the second column normalized by the average concentration at the corresponding axial positions, C cup (z). The fourth column shows the streamlines of the transverse flows with the reactive interface shown in green. See movie 3, available from stacks.iop.org/NJP/11/075028/mmedia, for an animation of the evolution of the concentration distributions. Colourbar below plot shows concentration colour scale. onto the uniaxial curve is only visible for Pe = 10 2 . Again we expect Sh(Z) in the entrance region to fall on the uniaxial curve for all Pe. Plots at very short axial length show agreement with this supposition, although the data are very noisy at such small distances, and one must infer by the overlap of the scattered data that they exhibit the same behaviour. Assuming this behaviour, the scaled entrance length can be inferred from the graphs, and is again a function of Pe, which decreases with increasing Pe. The asymptotic value of the Sherwood number, Sh ∞ trans , again increases with increasing Pe. As compared with figure 6(a), however, the inferred L ent /PeH decreases and Sh ∞ trans increases significantly more quickly with increasing Pe than in the case of transport to the stationary interface. This stronger dependence on Pe leads to very large values for Sh ∞ trans at high Pe, nearly 20 times larger than Sh ∞ stat and over 200 times larger than the asymptotic value in the absence of transverse flow.
We note several peculiarities in the curves in figure 8(a). All of the curves show some local variations in Sh(Z) that deserve comment. We attribute the non-periodic variations early in the entrance region and around the asymptotic value to the finite number of reaction events in these regions: the smoothness of the evaluation of Sh(Z) as in (3) is limited by the total number of particles and by the number of points at which Sh(Z) is evaluated. We also note coherent oscillations (e.g. Pe = 10 3 and 10 4 in figure 8(a)) and a significant rise in the average value of Sh(Z) for transport to the moving interface above the initial plateau (again for Pe = 10 3 and 10 4 in figure 8(a)). The spatial frequency of the oscillations equals the length of a halfcycle and arises due to an oscillation in the concentration of fluid that is swept down toward the reactive interface. The rise in Sh(Z) at large axial distance is attributable to the inverted boundary layer that forms in the case of a moving reactive interface, as discussed above in reference to figure 7: the persistent region of low concentration that forms opposite the reactive boundary significantly affects the average bulk concentration, C cup , used in calculating Sh(Z) as in (3). Yet, it is the difference between the concentration just outside the boundary layer and the concentration at the reactive interface itself that drives transfer; as the inverted boundary layer develops, this concentration becomes larger than C cup and leads to the rise in Sh(Z) to a second plateau. The value of Sh(Z) on the earlier plateau, before this rise in Sh(Z), is the value that will be tested in figure 11 against the theory presented in section 2.2.2, because the theoretical development lacks this detailed accounting for the actual shape of the concentration profile. This effect provides interesting insight for the design of flows for rapid transfer to a reactive interface: if a flow sends depleted fluid to a region far from the reactive interface, thereby maintaining a higher concentration just outside the boundary layer near the reactive interface relative to the bulk average concentration, transfer will be faster than in a flow with a similar surface process that immediately mixes its depleted fluid into the bulk (as in the case of the stationary boundary in figures 5 and 6). Indeed, the ideal flow would deplete the fluid near the reactive interface quickly and completely and then send that depleted fluid as far from the surface as possible, such that it only returns to the reactive interface once the entire bulk has passed near the surface and been depleted to a similar extent.
For the case of transport to the moving interface, there are again similarities and differences between the chaotic (figure 8(a)) and non-chaotic flows ( figure 8(b)). For the first few equivalent cycle lengths (z < 60H ), the curves in figure 8(b) are nearly identical to those in figure 8(a): Sh(Z) appears to reach the same asymptotic values in the non-chaotic flow as it did in the chaotic flow. But, for z > 60H , all curves collapse onto a new Pe-independent curve for Sh(Z) in the non-chaotic flow, which decays with increasing axial distance to an asymptotic value that is an order of magnitude larger than the uniaxial case. As in figure 6(b) for transport to the stationary interface, this decrease is due to the reinforcement of the boundary layer after depleted fluid has circled the vortices and returned to the reactive interface. Because the transfer to the moving interface is more efficient than that to the stationary interface, the depleted fluid that returns to the interface is significantly more depleted and results in a larger reinforcement. It is important to note, however, that the final plateau in Sh(Z) for the non-chaotic flow is still significantly larger than that for uniaxial flow. This increase is due to two factors: firstly, the flow near the moving interface is much more conducive to transport than the stationary boundary in the uniaxial case. Secondly, the final thickness of the boundary layer only extends to the cores of the vortices, rather than through the full depth of the channel; this distance is still independent of Pe, however, so we expect that all Pe will eventually collapse to this final plateau value. This collapse to a Pe-independent Sh(Z) again demonstrates a failure of the modified Graetz behaviour in the non-chaotic flow, as noted in reference to figure 6(b). The overall Sherwood number for transport to the moving interface averaged over the entire device will be much smaller than in the chaotic case, due to the fact that Sh(Z) decreases by an order of magnitude or more due to the reinforced depletion seen in figure 7.
The strong oscillations in Sh(Z) for the chaotic flow in figure 8(a) are absent for the nonchaotic flow in figure 8(b); this observation supports the idea that the oscillations are connected to the modulation of the chaotic flow. We again observe the rise in Sh(Z) for Pe = 10 3 at large axial distance, with a nearly identical shape to that in figure 8(a). We again attribute this increase to the formation of an inverted boundary layer (as seen in figure 7 and movie 2). This increase in Sh(Z) only appears at Pe = 10 3 because for Pe 10 4 , the asymptotic value of Sh(Z) approaches uniaxial; chaotically stirred with r = 1/3, u trans = U/5 and L 1/2 = 10H . The first column presents the concentration, C(x, y, z), in the cross section at a distance z downstream that is equivalent to twenty-five cycles of the mixer (z = 50L 1/2 = 500H ). The second column shows the width-averaged concentration in the cross section,C(y, z), as a function of axial distance, and the black line marks the position along the length of the channel corresponding to the crosssectional profile in the first column. The third column is equivalent to the second column normalized by the average concentration at the corresponding axial positions, C cup (z). The fourth column shows the streamlines of the transverse flow with the reactive interface shown in green. See movie 4, available from stacks.iop.org/NJP/11/075028/mmedia, for an animation of the evolution of the concentration distributions. Colourbar below plot shows concentration colour scale.
the Pe-independent curve for the non-chaotic flow, which has a lower asymptotic value than the corresponding value of Sh ∞ trans for the chaotic flow. Figure 9 shows the evolution of the concentration profile in a vertically symmetric flow, driven by a slipping boundary at the ceiling and floor as in figure 1(c). The uniaxial case looks similar to the analogous cases in figures 5 and 7. The entire device remains in the entrance region; the boundary layer continues to grow for the entire length of the channel. In the chaotic flow, the depletion of solute proceeds in two steps: firstly, transfer of solute to the bottom reactive interface rapidly depletes the lower half of the channel, and the concentration drops to a small, nearly constant value; secondly, once a significant gradient exists between the top and bottom Figure 10. Local Sherwood number, Sh(Z), as a function of scaled axial distance, Z = z/PeH , for transport to the green (moving) interface in the flow shown in figure 1(c) with r = 1/3, u trans = U/5 and L 1/2 = 10H . Colours signify Pe (10 2 (purple), 10 3 (red), 10 4 (green), 10 5 (blue), 10 6 (grey)) at constant u trans = U/5, and the black dots signify the collection of all Pe for u trans = 0 (uniaxial Poiseuille flow).

Transport to moving interfaces with simultaneous transport across diffusive internal interfaces
halves, transfer from top to bottom slowly depletes the top half. We note that the normalized side view demonstrates the convergence of the concentration profile in the chaotically stirred flow within the first ten cycles of the device, after which the concentration profile is self-similar and continues to decay in magnitude with increasing axial distance. Figure 10 shows the dependence of the local Sherwood number on scaled axial distance for transport to the moving interface in the chaotic flow shown in figure 9 and movie 4. As discussed above, transfer of solute to the bottom interface proceeds in two steps: firstly, the bottom half of the channel empties due to a rapid transfer process to the moving interface. This process defines the first plateau (at short axial distance) for each value of Pe in figure 10. This plateau corresponds to the asymptotic value of Sh(Z) predicted in table 1 for transport to the moving interface in a duct of square cross section. This step is similar to the transfer to the moving boundary in the channel of lower aspect ratio ( figure 1(a)), and, as such, the early behaviour of Sh(Z) shares some important features with figure 8(a), with slightly higher plateau values due to the difference in geometry. Secondly, after an axial distance, z ∼ = 40H , the top half of the channel empties slowly, as solute passes through both the diffusive internal interface defined by the plane of symmetry and the fluid-solid interface at the bottom, moving boundary. After this transition, Sh(Z) drops to a second, lower plateau; the two transfer steps in series-at the symmetry plane and at the fluid-solid boundary-define this asymptotic value.
The behaviour of Sh(Z) for transport to the reactive interface in a uniaxial flow with square cross section, shown in black in figure 10, demonstrates collapse for all Pe and an asymptotic value that is the same as that in the rectangular cross section with aspect ratio of W/H = 2.
The data for Sh(Z) for transport to the moving interface at Pe = 10 2 show an entrance region that falls on this curve as well, and so we conjecture that all values of Pe demonstrate the same entrance region behaviour at sufficiently short axial distances. The oscillations seen in figure 8(a) also appear in figure 10 for the chaotically stirred flow at Pe = 10 3 and 10 4 . The appearance of these oscillations makes sense as this flow is modulated, as explained for figure 8 above. The rise in Sh(Z) that is seen in the presence of both transverse flows in figures 8(a) and (b) is absent here. As the depleted fluid is swept slowly near the centreline, transport across the centreline reconcentrates this fluid, leading to a more uniform concentration in the lower half of the channel than that which is present in the flow driven by only the bottom boundary; the inverted boundary layer does not form.

Asymptotic values of Sherwood number and comparison to predicted correlations
Before moving on to compare the results for Sh ∞ from figures 6, 8 and 10 with the predictions in table 1, we must consider how the final asymptotic values in figure 10 are related to the Sherwood number at the centreline, Sh ∞ centre . The process that defines the final asymptotic value of Sh(Z) is a combination of Sh ∞ centre and Sh ∞ trans from table 1. The usual way to combine such processes is to apply the analogy of resistances in series, where the resistance to each transport process is the inverse of the Sherwood number for that process (Bird et al 1960). We note, however, that Sh ∞ centre is defined in terms of the difference in concentration across the centreline while Sh ∞ trans is defined in terms of the difference between the concentration in the lower half of the channel and the concentration at the reactive interface. The overall effective Sherwood number, Sh ∞ overall , is defined in terms of the overall average concentration in the entire cross section. By equating the expressions for flux at the reactive interface due to transport out of the lower half by Sh ∞ trans , due to transport out of the entire bulk by Sh ∞ overall , and due to transport from the upper half to the reactive interface across the equivalent resistance imposed by Sh ∞ centre and Sh ∞ trans , we have After some algebraic manipulation and inserting the predictions from table 1 for Sh ∞ centre and Sh ∞ trans , we find that and this Sh ∞ overall is the value that is measured by the second set of plateaus in figure 10. Figure 11 presents the asymptotic values of the Sherwood numbers, Sh ∞ extracted from the data in figures 6, 8 and 10. This summary allows us to demonstrate the relative magnitudes of Sh ∞ for the various cases, and to evaluate the success of the predicted correlations in table 1. All of the data show good agreement in scaling with the predicted correlations at high Pe trans . The data for Sh ∞ trans for transport to the moving interface in the flows in figures 1(a) and (c) also show excellent agreement in their prefactors (within 2%). The data for Sh ∞ trans in the flow in figure 1(c) begin to fall off of the predicted correlation for Pe trans < 10 4 , which was also predicted due to the influence of the shear-like component of the velocity at the moving boundary as the boundary layer grows thick enough to sample the shear flow (see section 2.2.2). The data for Sh ∞ stat and Sh ∞ overall exceed the predicted correlations by about 8 and 12%, respectively at high Pe trans . For Pe trans < 10 3 both Sh ∞ stat and Sh ∞ overall exceed their respective correlations by increasing amounts with decreasing Pe trans . The plateau in the uniaxial case (Sh ∞ ∼2.55) is a lower bound on the possible values of Sh(Z) for all flows in these geometries and therefore Sh ∞ stat and Sh ∞ overall cannot decrease beyond Sh ∞ ∼2.55 as Pe trans decreases.

Mixing length, return length, and the modified Graetz behaviour
We now present results for the mixing length and return length in the flows in figures 1(a) and (b) and discuss the implications of these results on the ability of flows to sustain high rates of mass transfer associated with the modified Graetz behaviour.
In figure 12, we test the scaling of the mixing length with Pe by scaling the axial distance of the curves of the axial velocity-weighted variance of concentration as defined in (23) for flows with an initial concentration difference between two streams of equal area in the cross section. The inlet concentration condition in each case consists of a stream of uniform high concentration occupying the upper half of the channel and a stream of zero concentration occupying the lower half of the channel; this initial distribution is representative of the distributions that are generated by reaction at the top or bottom interface. By scaling the axial distance by power laws in Pe with scaling exponents between 0 and 1, as well as by ln(Pe), we have found reasonable estimates of the dependence of mixing length on Pe. Figure 12(a) shows the result of scaling the axial distance in the chaotically stirred flow in figure 1(a) by Pe 1/6 . This scaling provides notably Figure 12. Decay of velocity weighted variance (23) with scaled axial distance for determination of scaling of mixing length with Pe: (a) chaotically stirred flow as in figure 1(a) with r = 1/3, u trans = U/5; L 1/2 = 10H with axial distance scaled by Pe 1/6 ; (b) non-chaotically stirred flow as in figure 1(b) with r = 1/2, u trans = U/5 with axial distance scaled by Pe. Colours signify Pe (10 3 (red), 10 4 (green), 10 5 (blue), 10 6 (grey)). Representative concentration profiles as a function of axial distance at Pe = 10 5 are shown to the right. better collapse than ln(Pe), that was expected due to Ranz (1979) and Villermaux et al (2008). Scaling by ln(Pe) causes the curves for various Pe to have similar shapes, i.e. similar curvatures and slopes, but the curves do not coincide or collapse to any significant extent. We conclude that this weak power law is a better representation of the scaling of the mixing length based on this velocity-weighted measure of mixing. The concentration profiles to the right in figure 12(a) show the mixing process at Pe = 10 5 .
For the non-chaotic flow in figure 1(b), the mixing curves in figure 12(b) demonstrate that mixing occurs in two steps that occur in series and present distinct scaling with respect to Pe. These two distinct processes can be seen qualitatively in the concentration profiles to the right of figure 12(b), where the transition between the two regimes occurs in the vicinity of z/H = 160.
The early process mixes the cores of the vortices and the exteriors of the vortices separately without significant mixing between these two regions. This process is analogous to the mixing of the depleted fluid from the boundary layer throughout the exterior of the vortices in the nonchaotic flow; this portion of the curves collapses when plotted against z/Pe 1/3 (characteristic of transfer across a shear flow). The second process is responsible for mixing the fluid in the cores of the vortices with the fluid outside of these cores. This second process is much slower than the initial process, with a mixing length that scales linearly with Pe, as demonstrated in figure 12(b). It is this second process on which we depend for mixing of the high concentration fluid in the cores of the vortices with the depleted fluid that has been swept off of the reactive interface; therefore the scaling of the relevant mixing process for transport to reactive interfaces in the non-chaotic flow is linear in Pe.
To determine the scaling of return length with Pe for a given flow, we injected a tracer particle into the flow and recorded the axial positions of its first 100 000 collisions with each reactive interface. Return length is defined as the axial distance traversed between subsequent collisions (figure 3). The distributions of return lengths were calculated and plotted, as shown in figures 13(a) and (b). The distribution of return lengths for returning to the stationary interface in the chaotic and non-chaotic flows are continuous, relatively smooth and both have the same general features, as seen at Pe = 10 4 in figure 13(a): an initial decay of short distance return lengths (this segment of the histogram contains most of the calculated trajectories), a minimum at intermediate return length, and a wide peak of long return lengths. The initial decays are nearly identical for both cases; the minima occur at similar return length, and the peak of long return lengths are of essentially the same shape. The initial decay corresponds to the growth of the boundary layer, wherein a particle collides with the reactive interface many times while remaining near the interface. This behaviour is seen in trajectories of particles that have clusters of short return lengths; they appear to skip across the interface like a rock skipping across the surface of a lake. When the particle is pulled away from the reactive interface into the bulk of the flow, there is a decrease in the probability of encountering the interface until the flow brings the particle back around to the interface again. This decrease in the probability of collision corresponds to the minimum in the distribution of return lengths. Finally, the wide peak after this minimum corresponds to the distribution of return lengths that involve a trip into the bulk prior to returning. We are interested in the ability of the flow to mix fluid that has been removed from the reactive interface before it returns from the bulk to the interface, and accordingly we only use the return lengths beyond the minimum in the distribution to calculate the characteristic return lengths that are plotted in figure 13(c). Figure 13(b) shows the distribution of return lengths for a particle returning to the green, (moving) reactive interface in figures 1(a) and (b). Again the distributions for particles in the chaotic and non-chaotic flows are quite similar to each other, although the features of the distributions are somewhat different from those in figure 13(a). The distributions have an initial decay that corresponds to the growth of the boundary layer as in figure 13(a). This portion of the distribution contains far fewer of the return lengths than the initial decay for return to the stationary surface. The time spent traversing the interface is much smaller, leaving less time during which to collide with the surface before being swept into the bulk. The minima in the distributions of return lengths, in this case, become gaps; no trajectories in our sample had return lengths between the end of the initial decay and the beginning of the final peak. While it is possible to spend multiple cycles near the stationary interface without being swept definitively into the bulk, the transverse flow always sweeps particles near the moving interface away within with Pe = 10 4 : (blue curve) r = 1/2, u trans = U/5; (grey curve) r = 1/3, u trans = U/5, L 1/2 = 10H . (c) Dependence of characteristic return length on Pe for return to the reactive interface: (crosses) r = 1/3, u trans = U/5, L 1/2 = 10H ; (circles) r = 1/2, u trans = U/5. Colours: (red symbols) return to the blue (stationary) interface; (green symbols) return to the green (moving) interface; solid line = 6Pe 2/3 ; dashed line = 4.5Pe 1/2 . one half-cycle of the mixer due to the magnitude of velocities and the length of the half cycle. The distributions again show a wide peak at large return length, and we use this portion of the distribution for calculation of the characteristic lengths plotted in figure 13(c) as stated above.
The characteristic return lengths plotted in figure 13(c) demonstrate scaling of Pe 2/3 for return to the stationary interface and Pe 1/2 for return to the moving interface. These results validate the simple arguments in section 2.3.2 and in (24) based on the behaviour of the average concentration at long axial distance and its relation to the Sherwood number for transport to the reactive interface. The chaotic character of the flow has no effect on the scaling of the return length.
We can now consider the ability of the flows considered here to maintain high rates of mass transfer in the framework of the modified Graetz behaviour out to arbitrary axial lengths and values of Pe. The chaotic flow in figure 1(a) presents modified Graetz behaviour in the ranges of Pe and of axial distances considered in this work for transfer to both stationary and moving interfaces (figures 6(a) and 8(a)). Based on the discussion in section 2.3, this observation implies that, at each Pe considered, the mixing length is shorter than the return length, such that the depleted fluid is mixed with the bulk before returning to the reactive interface. Thus, for all Pe considered in this work, we predict that the chaotic flow will behave according to the modified Graetz behaviour out to arbitrarily large axial distances. Having shown that the return length scales more strongly than the mixing length with increasing Pe for both reactive interfaces (stationary and moving-figure 13(c)), we also induce that the modified Graetz behaviour will hold out to arbitrarily large Pe. The non-chaotic flow in figure 1(b) demonstrates modified Graetz behaviour for transport to the stationary surface for Pe < 10 5 over at least two orders of magnitude in axial distance (figure 6(b)); this observation leads us to conclude that for low to intermediate Pe, the non-chaotic flow is able to transport solute out of the cores of its vortices fast enough to keep up with transport to the stationary interface, and will therefore show modified Graetz behaviour out to arbitrarily large axial distance. Yet, for transport to the stationary interface at large Pe (figure 6(b)) and for transport to the moving interface at all Pe considered ( figure 8(b)), the modified Graetz behaviour breaks down at sufficiently large axial distance. Having shown that the relevant mixing length scales more strongly than the return length with increasing Pe for both reactive interfaces, we induce that the modified Graetz behaviour will not hold for the non-chaotic flow for any Pe larger than those considered in this work.

Further generalization of the modified Graetz behaviour
The transport processes involving the convectively disconnected sets in the chaotic flow of figure 1(c) and figures 9 and 10 are a special case of a more general property of the modified Graetz behaviour. In the case of a symmetric flow as in figure 1(c), the internal interface presents a well-defined, purely diffusive boundary to transport between the two disconnected sets. Generally speaking, the presence of any poorly mixed islands in the bulk or on the wall can have similar effects on the overall transport process. If each set is well-mixed individually, all sets should show modified Graetz behaviour and a cascade of concentration as seen in figures 9 and 10. Any sets that are not well mixed individually and any interfaces that lead to a continuously broadening boundary layer for transport across the interface will show standard Graetz behaviour with an asymptotic value, Sh ∞ , that is determined primarily by the absolute size of the set or the maximum thickness of the boundary layer, as in the behaviour of the non-chaotic flow at large axial distance (see figures 5-8). The relative size of the sets and the properties of the velocity conditions near the interfaces between sets will determine how strongly each set affects the overall Sherwood number.
As an example of the generality of this behaviour figure 14 shows Sh(Z) and representative concentration profiles for a flow similar to the flow in figure 1(a), but with a width equal to the height, W = H , rather than W = 2H in figure 1(a). We see in the concentration profiles Figure 14. Local Sherwood number, Sh(Z), as a function of scaled axial distance, Z = z/PeH , for transport to the green (moving) interface in the flow shown in figure 1(a) with r = 1/3, u trans = U/5, and L 1/2 = 10H , but with W = H . Colours signify Pe (10 2 (purple), 10 3 (red), 10 4 (green), 10 5 (blue), 10 6 (grey)) at constant u trans = U/5, and the black dots signify the collection of all Pe for u trans = 0 (uniaxial Poiseuille flow). Inset: streamlines of the transverse flow, with reactive interface in green. Right: concentration profiles at several axial distances for each Pe in the plot of Sh(Z). that a fraction of the cross section is well mixed, but there is also a region near the top wall that is poorly mixed with the rest of the bulk. As solute reacts out of the system through the green (moving) interface, the average concentration in the well-mixed set decreases quickly while the poorly mixed set remains at roughly the inlet concentration. As in figure 10, this process corresponds to the initial plateaus seen in Sh(Z) in figure 14. At large axial distances, transport from the poorly mixed set near the top wall across the well-mixed set to the bottom (moving) interface defines the final asymptotic value of Sh(Z). It is interesting to note that the spacing between the initial and final plateau values of Sh(Z) in figure 10 was independent of Pe, implying a constant multiplicative factor between Sh ∞ overall and Sh ∞ trans and the same scaling with Pe for Sh ∞ centre and Sh ∞ trans , as predicted in table 1 and shown in figure 11. In figure 14, the spacing between the plateau values at Pe = 10 3 appears to be larger than that at Pe = 10 4 , and the second plateau is not even seen at Pe = 10 5 (although it may exist at larger axial distance), implying a scaling for Sh ∞ overall (second plateau) that is different from (in fact higher than) the scaling for Sh ∞ trans (first plateau). A higher scaling for Sh ∞ overall also implies a higher scaling for the Sherwood number for transfer between the two disconnected sets. This increased scaling makes sense, as there is some convective transfer across the interface between the sets in figure 14 at the end of each half-cycle, while there is only diffusive transfer across the well-defined interface at the centreline in figures 9 and 10.

Conclusions
In this work, we have extended the concept of the modified Graetz behaviour to the problems of mass transfer to moving interfaces and internal diffusive interfaces. Our predicted correlations (see table 1) for the scaling and prefactor of the asymptotic value of the Sherwood number for mass transfer to stationary, moving, and diffusive interfaces accurately represent the observed results from numerical simulations (see figure 11). Concentration profiles in the case of transport to moving interfaces revealed an interesting property of the LDC flows with which we work: as the depleted fluid is swept off of the moving reactive interface, in some cases it forms an inverted boundary layer far from the reactive interface. This inverted boundary layer sequesters depleted fluid to a large region near the stationary boundary opposite the moving, reactive interface; this segregation maintains high concentration near the reactive interface and allows for higher rates of mass transfer than would be achieved were the concentration in the bulk uniform. Although non-chaotic, 3D flows initially show higher Sh(Z ) relative to uniaxial flow, they eventually lead to an asymptotic value, Sh ∞ , that is independent of Pe, as in the standard Graetz case.
We have also presented the criterion for the existence of modified Graetz behaviour: for a flow at a given Pe, the flow must mix depleted fluid with the bulk before it returns to the reactive interface; the characteristic length for return of a depleted fluid element (a 'depleton') to the interface must be longer than the length required to homogenize it into the bulk. We argued further that the persistence of modified Graetz with increasing axial length and Pe should depend on the Pe-scaling of the return and mixing lengths: if the return length grows more rapidly than the mixing length, then, if modified Graetz is observed at finite Pe and length, it should persist for all higher values of these parameters. By simulation, we found weak scaling of the mixing length (∼Pe 1/6 ) relative to return length (∼Pe 2/3 for return to the stationary interface and ∼Pe 1/2 for return to the moving interface) in the chaotically stirred flow. We thus conclude that modified Graetz should persist in this flow; this conclusion is borne out by our direct calculation of Sh from simulations out to Pe = 10 6 . In contrast, for the non-chaotic flow, strong scaling of the mixing length (∼Pe) should lead to the breakdown of modified Graetz at high Pe; this prediction is consistent with our observation that non-chaotic flow fails to sustain the modified Graetz behaviour beyond intermediate Pe.
Finally, we note that there is not a direct analogy between the results for transfer to moving interfaces in the LDC flows and for transfer to the stationary grooved surfaces that drive similar flows in the SHM. As such, we cannot expect to see the same increase in Sh ∞ in the SHM that we report here for Sh ∞ at the moving interface. We are currently pursuing further investigations of transport to grooved boundaries.