Shear-flow dispersion in turbulent jets

We investigate the transport of a passive scalar in a fully developed turbulent axisymmetric jet at a Reynolds number of $\mathit{Re}=4815$ using data from direct numerical simulation. In particular, we simulate the response of the concentration field to an instantaneous variation of the scalar flux at the source. To analyse the time evolution of this statistically unsteady process we take an ensemble average over 16 independent simulations. We find that the evolution of $C_{m}(z,t)$ , the radial integral of the ensemble-averaged concentration, is a self-similar process, with the front position and spread both scaling as $\sqrt{t}$ . The longitudinal mixing of $C_{m}$ is shown to be primarily caused by shear-flow dispersion. Using the approach developed by Craske & van Reeuwijk (J. Fluid Mech., vol. 763, 2014, pp. 538–566), the classical theory for shear-flow dispersion is applied to turbulent jets to obtain a closure that couples the integral scalar flux to the integral concentration $C_{m}$ . Model predictions using the dispersion closure are in good agreement with the simulation data. Application of the dispersion closure to a two-dimensional jet results in an integral transport equation that is fully consistent with that of Landel et al. (J. Fluid Mech., vol. 711, 2012, pp. 212–258).


Introduction
Shear-flow dispersion (Taylor 1953(Taylor , 1954bAris 1956) is one of the primary sources of mixing in integral models for passive scalar transport and as such has wide-ranging practical applications, including contaminant transport in the atmosphere and ocean, nutrient delivery, the spread of smoke from fires and the discharge of waste effluent in streams (see, e.g. Fischer et al. 1979). Caused by lateral (cross-stream) gradients in a mean velocity, shear-flow dispersion was first identified by Taylor (1953), who examined the transport of a solute in both laminar (Taylor 1953) and turbulent (Taylor 1954b) pipe flow. It was demonstrated that for sufficiently large times (Taylor 1954a) the effective longitudinal (streamwise) mixing is determined by a balance between longitudinal advection and radial mixing.
The canonical example of shear-flow dispersion is the release of a passive scalar into a bounded one-dimensional flow in a pipe or a plane channel (see figure 1). At t = 0, one assumes that the mean scalar concentration c, where the overline denotes an ensemble average, is uniformly distributed over the pipe for z 0 and equal to zero for z > 0. One is typically interested in the integral concentration where the upper limit of integration r d corresponds to the radial extent of the flow. Hence, for a pipe of radius r m , r d = r m and C m = r 2 m c , where denotes an average over the cross-section of the pipe. At t = 0 the dependence of C m on z corresponds to a Heaviside step function. For time t > 0 the step change in C m will become progressively smoother due to three distinct, yet closely related, physical processes. D1 The diffusive flux −κ∂ z c (molecular diffusion). D2 The ensemble covariance wc − w c (turbulent mixing). D3 The spatial covariance w c − w c (shear-flow dispersion).

Shear-flow dispersion in turbulent jets
Here w is the longitudinal velocity and κ is the molecular diffusion coefficient. Molecular diffusion is a manifestation of an averaging operation over individual particles in the fluid and is expected to play a negligible role at large scales in high-Reynolds-number flows. Turbulent mixing is a result of the chaotic motion of the flow, the covariance wc − w c corresponding to a turbulent scalar flux. Shear-flow dispersion, on the other hand, is the result of a correlation of the average velocity w with the average concentration c over space, and enters the problem because the integral concentration C m is the primary unknown. When regions of high mean concentration coincide spatially with regions of high longitudinal mean velocity (see figure 1), w c = w c , and one finds a dispersive flux. It should be noticed that the turbulent scalar flux and the dispersive flux in D2 and D3 are defined relative to the ensemble average scalar flux, w c, and the spatially averaged flux, w c , respectively. In many large-scale practical applications, such as open channel flow, one finds D3 D2 D1 (see e.g. Elder 1959). Moreover, Taylor (1953Taylor ( , 1954b) demonstrated that longitudinal dispersion is inversely proportional to the turbulent diffusion coefficient associated with the flow, indicating that mixing by turbulence inhibits shear-flow dispersion. Further details can be found in appendix A, in which we summarise the dispersion theory of Taylor (1953) in the context of pipe flow. For a general introduction to dispersion theory the reader is referred to Chatwin & Allen (1985).
Taylor's original analysis of pipe-flow dispersion was subsequently generalised to pipes of arbitrary cross-section by Aris (1956), who obtained solutions for the longitudinal moments of the solute concentration and demonstrated that the longitudinal dispersion coefficient and longitudinal turbulent diffusion coefficient are additive. Gill (1967) later considered a general series expansion, involving higher longitudinal derivatives of the mean concentration, to determine the radial dependence of the concentration. Subsequent work has addressed the way in which the concentration profile approaches a Gaussian form (Chatwin 1970), the asymptotic behaviour for small times (Chatwin 1977) and the concentration at large distances from its centre of mass (Haynes & Vanneste 2014). Following the work of Brenner (1980b), a generalised Taylor dispersion theory emerged, which was not restricted to the unidirectional flows on which previous studies had focused. Generalised Taylor dispersion theory has been applied in a wide range of fields, including sedimentation (Brenner 1979), flows through porous media (Brenner 1980a) and chemically reacting flows (Shapiro & Brenner 1986). However, in spite of the many refinements that have been made to Taylor's theory of dispersion it has not, to the best of the authors' knowledge, been applied explicitly in the analysis of scalar dispersion in jets. Very recently, dispersion theory was used to develop a closure for longitudinal mixing in statistically unsteady jets (Craske & van Reeuwijk 2015b). There, the closure was applied to the mean kinetic energy flux and produced results that were in good agreement with direct numerical simulation. In this work we study the canonical problem of scalar dispersion in statistically steady jets and show that the dispersion theory developed in Craske & van Reeuwijk (2015b) provides an accurate representation of the dominant mixing processes.
Passive scalar transport in turbulent jets has received significant attention in the literature (e.g. Paranthoen et al. 1988;Tong & Warhaft 1995;Warhaft 2000). However, the primary focus has been on steady releases; unsteady releases, for which shear-flow dispersion is expected to be significant, have received relatively little attention. An exception is the work by Landel, Caulfield & Woods (2012), who investigated unsteady scalar transport in a two-dimensional jet, created in the laboratory by confining the jet to a narrow gap between two parallel plates. They observed that the core region of the jet, which is primarily associated with longitudinal advective transport, is surrounded by eddies that are responsible for mixing. Moreover, this structure was observed to be self-similar with height, which means that when normalised by suitable length and velocity scales the behaviour of the jet is independent of its streamwise coordinate. Longitudinal mixing was attributed to the stretching of fluid that occurs between the core and the eddies. Using a mixing-length assumption, justified by the self-similarity of the flow, Landel et al. (2012) formulated an advection-diffusion equation for the transport of the integral C m (z, t) of a scalar given by Here, z is the longitudinal (streamwise) coordinate, M m is the jet momentum flux, and K a and K d are empirical dimensionless parameters accounting for advection and eddy diffusivity respectively. z r FIGURE 2. Definition sketch. The grey area represents regions to which the scalar has spread at a particular time t and w is the ensemble average longitudinal velocity.
In spite of the fact that mixing-length models can provide only a limited representation of turbulence (Pope 2000), Landel et al. (2012) show that (1.2) provides an accurate means of predicting integral scalar concentrations in two-dimensional jets. However, it is neither clear what physical process determines the value of K a and K d , nor possible to state a priori which parts of the model can be applied to axisymmetric jets. Furthermore, it is unclear whether the mixing is caused by the longitudinal turbulent scalar flux or shear-flow dispersion (or both). In this work we show that for axisymmetric jets, the longitudinal turbulent scalar flux is not able to account for the observed mixing of the scalar integral, and that it is shear-flow dispersion that is the predominant cause of the mixing.
The paper is structured as follows. In § 2 we define the problem, introducing the governing equations and the variables that are of primary interest. In § 3 we describe the direct numerical simulation (hereafter DNS) used to obtain the data and describe how the ensemble average was constructed. In § 4 we demonstrate that the scalar transport can be viewed as a self-similar process by analysing the propagation speed of the scalar and its longitudinal spreading rate. We explain how Taylor's model for shear-flow dispersion in pipes (Taylor 1953) can be applied to jets in § 5, and compare the resulting dispersion closure with DNS observations in § 6. We apply our dispersion closure to planar jets in § 7 (additional details are provided in appendix B), which reveals that it is in one-to-one agreement with the model proposed by Landel et al. (2012). A summary of the findings and ideas for future work are provided in § 8.

Problem definition
A schematic for the flow considered is provided in figure 2. The turbulent jet is statistically steady, axisymmetric and produces a mean velocity field (u, w) in coordinates (r, z), where r denotes the radial and z the longitudinal coordinate. For t < 0, the scalar concentration is equal to zero everywhere. For t 0, a continuous release of scalar is activated at the inlet (hereafter referred to as the source), which produces a front that propagates and spreads as it travels along the jet. Due to the non-uniform velocity distribution, the scalar will be advected at different rates depending on its position in the jet: fluid elements in the core regions will tend to travel faster than elements in the outer region, thereby modifying the radial distribution of the scalar downstream. The motion of the jet is governed by the incompressible Navier-Stokes equations: We use cylindrical coordinates (r, φ, z), where r is a radial coordinate, φ is an azimuthal coordinate and z is a longitudinal coordinate, with corresponding velocity field u(r, φ, z, t) ≡ (u, v, w). The fluid is of constant uniform density ρ, the pressure relative to a hydrostatic balance is denoted by p and the kinematic viscosity of the fluid is ν. The passive scalar is governed by the transport equation where κ is the molecular diffusion coefficient. In this work the velocity field is statistically stationary and, due to the fact that the mean azimuthal component of velocity is equal to zero, swirl free. It is therefore convenient to define the time and azimuthal average velocity u ≡ (u(r, z), 0, w(r, z)). For the scalar concentration c, which is statistically unsteady, we restrict ourselves to an ensemble and azimuthal average c(r, z, t).
The focus of the present study is on the behaviour of integrals over lateral slices of the jet. In particular, we define the volume flux Q m (z) and the mean momentum flux M m (z): respectively, where r d denotes the radial extent of the jet, which will be defined precisely in § 3. In addition to the integral scalar concentration C m defined in (1.1), we define the integral mean scalar flux F m and the integral turbulent scalar flux F f , according to Here, χ denotes an ensemble average and χ a fluctuation from the ensemble average such that χ = 0. Formally, the integral quantities above are defined in the limit r d → ∞, but from a practical perspective it is necessary to integrate to a finite limit. This avoids complications associated with integration limits at infinity (see e.g. Kotsovinos 1978), allowing one to focus on the fluxes in the jet, rather than the induced ambient flow. The volume flux Q m and momentum flux M m can be used to define characteristic length and velocity scales r m ≡ Q m /M 1/2 m and w m ≡ M m /Q m respectively. We define a mean concentration c m for the flow from the integral concentration C m according to c m ≡ C m /r 2 m . Taking a time and ensemble average of (2.1) and (2.2) for Re ≡ 2M 1/2 m /ν 1, and integrating over lateral slices of the jet yields the integral jet equations (see e.g. which express volume conservation and momentum conservation respectively, and α is the classical entrainment coefficient (see e.g. Morton, Taylor & Turner 1956).
It is immediately apparent from (2.6a,b) that Q m = 2αM 1/2 m (z − z v ) for a conserved momentum flux M m and a virtual origin z v , the latter corresponding to the location at which the volume flux is equal to zero.
When molecular transport in the longitudinal direction is neglected, the ensemble average of (2.3) can be integrated over lateral slices of the jet to give which is an integral conservation equation for the scalar and is the principal focus of this work. In particular, our discussion of the problem will focus on the mean and turbulent dimensionless scalar flux, defined as respectively.

Simulation details
We simulate an axisymmetric turbulent jet driven by an isolated source of steady momentum flux M 0 and volume flux Q 0 . The source is located at the centre of the base of a cuboidal domain of size L x L y L z = 44 2 × 66 characteristic source radii, r 0 , where r 0 ≡ Q 0 /M 1/2 0 . The motion of the jet is governed by the incompressible Navier-Stokes equations (2.1) and (2.2) in a statistically steady state. To release a passive scalar into the flow at t = 0 we impose a step change in the scalar flux at the source: where H is the Heaviside step function. Equations (2.1)-(2.3) are solved numerically using N x N y N z = 768 2 × 1152 computational cells over a Cartesian grid. The code for the DNS of (2.1)-(2.3) employs a spatial discretisation of fourth-order accuracy that conserves volume, momentum and energy. Integration in time is performed using a third-order Adams-Bashforth scheme (further details can be found in Craske & van Reeuwijk 2015a). On the vertical and top faces of the domain we impose open boundary conditions that allow fluid to enter and leave the domain without modifying the entrainment rate of the jet (Craske & van Reeuwijk 2013). The approximately circular source at the base of the domain is obstructed by two narrow orthogonal strips of width ≈2r 0 /5 whose effect is to reduce the length of the potential core (see e.g. Craske & van Reeuwijk 2015a). We initiate the turbulence by applying uncorrelated perturbations of 1 % to the velocities in the first cell above the source. A constant mean scalar flux is maintained by imposing a Dirichlet boundary condition on both w and c at the source. Data for a fully developed statistically steady jet were acquired by running simulations for a total duration of not less than 6000τ 0 , where τ 0 ≡ Q 2 0 /M 3/2 0 is a characteristic time scale relevant to the source. Data from the first 3000τ 0 time units were discarded to eliminate large-scale transient behaviour, before statistics were obtained over the remaining duration of 3000τ 0 . To confirm that the integral scales of the turbulent jet are not influenced by viscous effects we simulate two jets, L and H, with Reynolds numbers Re ≡ 2M 1/2 0 /ν of 4815 and 6810 respectively. For each of these cases we impose a constant flux F 0 of scalar at the source in order 34 J. Craske, A. L. R. Debugne and M. van Reeuwijk to determine the steady-state concentration statistics. To ensure that the results are independent of grid resolution we perform a simulation H * , of reduced resolution, that uses 512 2 × 768 computational cells, but in other respects is identical to H. In practice, dimensional parameters such as ν and κ are defined by specifying the dimensionless parameters Re and Sc, where Sc is the Schmidt number. Consequently, it is convenient to let the inlet velocity be equal to 1 and specify ν = 2M 1/2 0 /Re and κ = ν/Sc. Further details of the simulations can be found in table 1.
Azimuthally averaged data were obtained by partitioning the domain into concentric cylindrical cells and averaging over all cells lying within a given shell. Therefore, in the steady-state simulations the mean χ of a quantity χ is obtained from an azimuthal average, in addition to a time average, and is consequently a function of r and z. To compute integrals such as C m , defined in (1.1), and the fluxes Q m , M m , F m and F f (see § 2), we define the upper limit of integration r d according to w(r d , z, t) = 0.01 w m (z, t) (a similar criterion was employed in the analysis of experimental data by Kotsovinos & List 1977).
During the course of the steady simulations we write a set of 16 complete threedimensional field files to disk at time intervals of approximately 250τ 0 . These field files are used to obtain independent initial velocity conditions for the simulation of 16 statistically unsteady scalar transport simulations, which we shall refer to as U1-16. For the scalar field we impose an initial condition of zero concentration uniformly throughout the domain. An unsteady scalar transport simulation is then created by imposing a constant flux F 0 of scalar at the source.
During each individual unsteady simulation, rather than outputting discretely sampled instantaneous data, a time average is taken over an interval in time that is much smaller than the characteristic time scale of the problem. Specifically, we employ a time interval not exceeding 2τ 0 , which prevents artificial smoothing of the front that time averaging would otherwise cause. In the unsteady simulations the azimuthal direction remains statistically homogeneous, which enables us to take azimuthal averages. In addition, the reliability of the statistics was further improved by performing an ensemble average over the 16 unsteady simulations. Therefore, a mean quantity χ obtained from the unsteady simulations is a function of r, z and t.
The dimensionless radial profiles of the steady scalar field c and the turbulence fluxes u c and w c are shown in figure 3. The profiles were taken from the range z/r 0 ∈ [28, 55] and are convincingly self-similar. The comparison between cases L and H in figure 3 supports the view that the large-scale integral statistics are not significantly influenced by the Reynolds number, and their comparison to case H * demonstrates that the results are independent of the spatial grid resolution. For a validation of the velocity data, including Reynolds stresses and a comparison of Q m and M m to plume theory, the reader is referred to Craske & van Reeuwijk (2015a). From the steady-state simulation data we infer that the entrainment rate Having demonstrated that the integral statistics such as w c are approximately independent of the Reynolds number Re, we choose to use case L at Re = 4815, rather than case H with Re = 6810, to obtain the unsteady data, which ensures that the scalar field has ample resolution and that the maximum cell Péclet number, max(u · x)/κ, for grid spacing x, remains small.

Propagation and spreading of the scalar front
4.1. Scaling Following a step change in the source scalar flux, a disturbance in the integral concentration C m propagates and spreads in the longitudinal direction in the form of a front. The objective of this section is to track the front position, which we denote z * (t), and determine the rate at which it spreads. Figure 4 displays the scalar field at several instants in time. On the left-hand side of each panel is a vertical slice through the scalar field of a single member of the ensemble and on the right-hand side is the ensemble and azimuthally averaged scalar field. As the scalar is introduced into the domain it is mixed by the turbulent jet and transported in the positive z direction. At each time the scalar field has a similar form, in which the widest part of a given isoregion lies approximately midway between the origin and its leading edge in z.
A robust criterion for identifying the location of the scalar front is to define z * implicitly as where F * m ≡ F m0 /2 is the arithmetic mean of the mean steady-state scalar flux F m0 before and after the step change. The value of F m0 was obtained from steady-state data by averaging F m (z) over the interval z/r 0 ∈ [28, 55], and is slightly lower than F 0 due to the contribution from the turbulent scalar flux (i.e. F 0 = F m + F f ). In addition, we define the positions z 1/4 (t) and z 3/4 (t), which correspond to the longitudinal distances at which F m = F m0 /4 and F m = 3F m0 /4 respectively, and will allow us to quantify the longitudinal spread of the scalar. Figure 5 shows typical profiles of C m and F m obtained from the ensemble-averaged simulation data, with circles marking the location of the front z * . In the region r 0 z z * the concentration field is in a quasi-steady state and the integral concentration C m increases linearly in z to a maximum behind the front (figure 5a). Ahead of its Left-hand side: instantaneous radial slice; right-hand side: azimuthal and ensemble average. maximum, C m tends to zero smoothly with increasing z/z * . It is observed that the propagation speed of the front decreases with increasing t. Evident in the behaviour of both C m and F m is that the longitudinal extent of the front increases in time, which indicates the presence of some form of longitudinal mixing. Due to the fact that the velocity in the jet w m ∼ 1/z, one expects disturbances to propagate according to z * ∼ √ t. To investigate this we plot (z * − z v ) 2 against t − t v in figure 6, where t v can be regarded as the location of the virtual source in time. Using the steady-state value of z v , t v = −1.5τ 0 was obtained by fitting a straight line to (z − z v ) 2 , whose agreement with the simulation data in figure 6 confirms the expected scaling relation. Both z v /r 0 and t v /τ 0 are small in comparison to the size of the domain and the total duration of the simulations respectively, and consequently play an insignificant role in the far-field scaling observed in figure 6. Interestingly, z * ∼ √ t corresponds to the classical dispersion scaling, although for planar jets and plumes, for example, in which w m ∼1/z, one would expect to find a different power-law scaling. It is noteworthy that the propagation rate of the front appears to be slightly greater than w m (figure 6, dashed line). Furthermore, from figure 6, z 1/4 and z 3/4 can be seen to scale in proportion to √ t, which suggests that the scaling associated with the position of the front is identical to that associated with its longitudinal rate of spread, and that the process is therefore self-similar. The scalar is released continuously at the source with a constant source flux F 0 , therefore ∞ 0 C m dz = F 0 t. (4.2) Since C m ∝ z in the steady state, we conclude that we can normalise the profiles of C m by F * m t/z * in order to observe self-similarity. Profiles of C m normalised in this way, plotted against the normalised streamwise coordinate z/z * , are shown in figure 7(a). The approximate collapse indicates that, following a step change in the flux of concentration at the source, the transport of a passive scalar in a jet is indeed a self-similar process that can be described by the variable z/z * .

Longitudinal mixing
To quantify the mixing at the front, we fit the data in figure 5 In an effort to account for the observed level of mixing, we assume that the turbulent Schmidt number is approximately unity and compare the value of D e with the uniform turbulent viscosity ν T obtained from a gradient-diffusion hypothesis. Based on the measurements of Hussein, Capp & George (1994), Pope (2000) reports that in jets ν T ≈ 0.033 M 1/2 m . (More precisely, Pope (2000) expresses ν T in terms of the centreline velocity w c and the radial location r 1/2 at which w = w c /2 (the jet half-width) as ν T ≈ 0.028 w c r 1/2 . In a Gaussian jet w c = 2w m and r 1/2 = √ ln 2/2 r m , hence ν T ≈ 0.033 w m r m .) The significant difference between this estimation and D e implies that turbulence alone cannot account for the mixing observed at the front. An intuitive way of understanding the difference between ν T and D e is that it implies that the longitudinal extent of the front (determined by D e ) exceeds the radius of the jet (determined by ν T ). We will demonstrate that shear-flow dispersion accounts for the difference by examining the way in which the concentration profile is distorted in the vicinity of the front. As described in § 1, longitudinal mixing of an integral quantity such as C m can occur due to a uniform increase in w c over the radius of the jet (turbulent mixing) or due to a distortion of the mean concentration profile c (shear-flow dispersion). Shear-flow dispersion results in a local increase in the correlation w c , involving only mean-flow quantities, and therefore in a local increase in the mean scalar flux. It is therefore appropriate to examine radial profiles of c and w in the vicinity of the front, which are shown alongside the fluxes w c and w c in figure 7. The profiles displayed in figure 7(b-d) are normalised using centreline values c c (z, t) ≡ c(0, z, t) and w c (z) ≡ w(0, z) of the mean concentration and the mean longitudinal velocity respectively. The central observation that can be made from figure 7(b-d) is that as z/z * increases the profiles become narrower. The scalar is transported by advection to regions of relatively large z/z * when it is located close to the centreline of the jet, where the longitudinal velocities are highest. Consequently, the scalar profile becomes more peaked with downstream distance, which is evident in figure 7(b). Figure 7(c,d) shows that the turbulent flux does not scale in exact proportion to the mean flux in the vicinity of the front. Rather, the turbulent flux appears to increase at large values of z/z * , although it continues to make a relatively small contribution to the total flux. To quantify precisely the relative contribution to the total scalar flux made by the mean and turbulent flow, we plot the dimensionless scalar fluxes θ m and θ f respectively with respect to the normalised downstream distance in figure 8. The grey curves displayed in figure 8 correspond to different times and exhibit an approximate collapse with respect to the variable z/z * , which is consistent with the self-similarity that was described in § 4.1. In the steady state the radial dependence of the velocity and concentration profiles is approximately independent of z, which implies that θ m is constant. Consequently, shear-flow dispersion is evident as a departure from the steady-state value of θ m and is seen in figure 8 to be significant over a large extent of the z/z * domain. Figure 8 also indicates that in a quasi-steady state, for which z/z * 1, the contribution of the dimensionless turbulent flux (θ f ) amounts to the expected 10 %, but rises to 20 % at the front, where z = z * . Figure 8 confirms that in the vicinity of z = z * the local distortion of the concentration profile c, evident in figure 7(b), acts in correlation with w to increase the dimensionless scalar flux θ m . In particular, the dimensionless scalar flux ranges from approximately 1.0 (for z/z * 1) to 2.0. It is also evident that θ f increases in the vicinity of z = z * . Indeed, since θ f is an integral quantity, in addition to accounting for a local increase in the magnitude of the turbulent flux w c (see figure 7d), it accounts for a narrowing of the spatial distribution of c . Hence, it is to be expected that the behaviours of θ m and θ f share similar features. We attribute the fact that θ f does not appear to remain in constant proportion to θ m to the relatively small increase in the magnitude of the turbulent flux w c . It should be noted that beyond z/z * , the fluxes w c, w c and the concentration C m used in their normalisation approach zero, which makes it difficult to obtain an accurate estimation of θ f and θ m from a finite dataset.

Shear-flow dispersion in turbulent jets
As pointed out previously, fluid located close to r = 0 is transported by relatively high longitudinal velocities, and can therefore propagate further than fluid located at the periphery of the jet. It is therefore natural that at large distances from the source one sees a narrowing of the scalar distribution (see figure 7b). At the leading edge of the scalar field, we can therefore expect the scalar profile to take the form of a spike at r = 0. For a Gaussian velocity profile, the largest velocity is found on the centreline and has a value of 2w m , which implies an upper bound of θ m = 2. Hence, the mean scalar flux F m cannot exceed 2C m w m . Furthermore, if it is assumed that dz * /dt ≈ w m (z * ) ∝ 1/z * , in accordance with figure 6, then fluid on the jet centreline will reach z/z * ≈ √ 2. At z/z * ≈ √ 2 in figure 8 we observe that the sum θ m + θ f ≈ 2. For z/z * > √ 2, we see a reduction in θ f and that θ m approaches a value of 2, as expected.

A model for shear-flow dispersion in jets
Although shear-flow dispersion was originally developed for pipe flow (Taylor 1953(Taylor , 1954b, Craske & van Reeuwijk (2015b) demonstrated that it could be usefully extended to boundary-free shear flow, such as a jet. In this section we will summarise the key features of the approach and obtain a model for the dispersion of a passive scalar. In contrast to a jet, the steady-state distribution of a scalar in a pipe is uniform and bounded in the radial direction. Like jets, however, both the velocity distribution and the scalar distribution are self-similar, which means that their radial dependence is invariant in the longitudinal direction. Fundamental to both dispersion in pipe flow and dispersion in jets is the idea that changes in the mean scalar concentration in the longitudinal direction can cause a local departure from self-similarity over the lateral (radial) dimension. Specifically, the departure from self-similarity in the scalar is caused by lateral gradients in the longitudinal velocity and, in correlation with the longitudinal velocity, results in a local increase in the scalar flux.
Taylor demonstrated that, to leading order, the radial dependence of the modification to the uniform concentration is independent of the longitudinal coordinate. Consequently, one can define a dimensionless function g 1 (η) to describe the shape of the perturbed concentration profile, where η ≡ r/r m . Using this notation, the basic dimensionless steady-state concentration will be denoted g 0 (η). The functions g 0 and g 1 apply equally well to the analysis of pipe flow and jet flow, although for pipe flow the steady-state scalar concentration is constant, and therefore g 0 does not depend on η. The aim is to establish an expression for the dimensionless concentration g(η) = g 0 + L 1 g 1 , (5.1) where L 1 determines the degree to which the concentration profile departs from a steady state, and is therefore expected to depend in some way on the behaviour of the scalar in the longitudinal direction. Evidently, if the dimensionless velocity profile w/w m = f 0 (η) is known, (5.1) can be multiplied by f 0 and integrated with respect to η 2 (i.e. over a lateral slice of the pipe flow or jet). The resulting integral corresponds to the dimensionless flux θ m ≡ F m /(w m C m ): We define θ 0 ≡ f 0 g 0 as the dimensionless steady-state scalar flux, and θ 1 ≡ f 0 g 1 as the additional flux arising from changes in the shape of g(η). Here, the operator is a suitable non-vanishing spatial average over η, which for jets we define in § 5.2. The ultimate aim of a dispersion closure is to obtain an expression for θ m in terms of the integral concentration C m and its longitudinal derivatives. Instead of working with the ensemble average scalar concentration c, Craske & van Reeuwijk (2015b) worked with the dimensionless concentration where c m0 (z) is the mean steady-state concentration and C m ≡ c m /c m0 . Consequently, in the steady state C is independent of z in both jets and pipes, making it possible to recast Taylor's approach to dispersion in pipe flow in a form that is applicable to jets. In the following sections we will substitute (5.3) into the longitudinal scalar transport equation to relate the unknown perturbed profile g 1 to ∂ z C m .

Shear-flow dispersion in a pipe
To begin, we describe the classical dispersion problem in a pipe flow using the notation described above. To understand the evolution of a step change in scalar concentration in a pipe of uniform radius r m , Taylor employed a frame of reference moving with the mean velocity w = w m . For pipes, for which θ 0 = 1, as the steady-state concentration is uniform, this coordinate transformation ensures that the equation governing the cross-sectional average concentration does not contain mean flux terms. More generally, it is necessary to move at the velocity associated with the steady-state scalar flux w m θ 0 to eliminate the mean flux terms. Taylor subtracted from the governing equation for scalar transport the evolution equation for the cross-sectional average concentration, and found, for large times, a leading-order balance between radial transport of the perturbed scalar profile and longitudinal advection (see appendix A for details): where κ T is the eddy diffusivity. In particular, (5.4) shows that the longitudinal gradient of the dimensionless average concentration C m provides the forcing for the perturbation in the concentration profile L 1 g 1 . The separable form of (5.4) allows one to infer that to within an arbitrary multiplicative constant, if it is assumed that κ T does not depend on η. Equations (5.4) and (5.5) demonstrate that the amplitude of the scalar perturbation is proportional to the longitudinal gradient of the dimensionless concentration, and that g 1 , and therefore θ 1 , can, in principle, be determined a priori by solving an ordinary differential equation when f 0 is known.

Shear-flow dispersion in an axisymmetric jet
To apply Taylor's dispersion theory to jets we assume that the decomposition C = C m (g 0 + L 1 g 1 ) provides a leading-order representation of the dimensionless concentration field. The use of the normalised radius η ≡ r/r m (z) allows us to map the conical geometry of the jet, for which r m = 2αz, onto a cylindrical geometry. Moreover, the use of the dimensionless concentration C ensures that in the steady state the amplitude of the scalar perturbation L 1 = 0, in spite of the fact that c m ∼ 1/z, and therefore allows us to focus on departures from this state. Due to the geometry of the jet, it is unlikely that the ordinary differential equation satisfied by the perturbation g 1 (η) is the same as that which is implied by (5.4) for pipes. However, since we are primarily interested in the integral quantity θ 1 ≡ f 0 g 1 , rather than the particular form of g 1 (η), we will not attempt to obtain an exact equation for g 1 (η) in this work.
In the absence of a well-defined edge, we define the cross-sectional average χ of a quantity χ in the jet according to such that w m = w . For axisymmetric jets, a constant eddy viscosity ν T can be expressed in terms of the entrainment coefficient according to ν T = αw m r m /3 (see Craske & van Reeuwijk 2015a). Assuming that the turbulent Schmidt number ν T /κ T is approximately unity, we therefore set κ T = αw m r m /3. Furthermore, the dimensionless mean scalar concentration C m can be defined in terms of integral quantities, such that C m ≡ C m /C m0 , where C m0 is an integral steady-state concentration. Consequently, substitution of κ T , C m and the linear dependence of the jet width on z, r m = 2αz, into (5.5) results in In the steady state, for which C m ∼ z, it is clear from (5.7) that L 1 = 0. Using (5.2), the dispersion closure for the dimensionless scalar flux in an unsteady jet can be expressed as As noted previously, instead of calculating g 0 and g 1 exactly, we treat θ 0 and θ 1 as model parameters whose value can be determined from observation. Notable, however, is the fact that θ 0 can be determined from the steady-state data.

Model prediction
Substitution of (5.8) into (2.7) results in an advection-dispersion equation given by Hereafter we will assume that z and t are coordinates relative to the location of a virtual source (z v , t v ). It is noteworthy that the advection parameter θ 0 + 6θ 1 (cf. K a in Landel et al. 2012) contains a contribution from both θ 0 and θ 1 , which are parameters that each have a particular physical significance and are precisely defined integrals. In particular, θ 0 depends on steady-state properties of the flow and can therefore be estimated a priori. Although θ m = θ 0 in the steady state, recognising that the longitudinal turbulent scalar flux θ f will affect the relevant advection velocity, we equate θ 0 with the observed steady-state value of θ m + θ f ≈ 1.0. To estimate the value of the parameter θ 1 , we compare the coefficient of the mixing term on the right-hand side of (6.1) with the mixing coefficient D e that was observed in § 4. Specifically, we set θ 1 = αD e /(3M 1/2 m ) = 0.017. The prediction of θ m that we obtain using (5.8) with (θ 0 , θ 1 ) = (1.0, 0.017) and the observed behaviour of C m is shown in figure 8.
Although the model appears to overpredict the magnitude of θ m , it shows a good agreement with the variation of θ m with respect to z/z * . In fact, the model shows a better agreement with the magnitude of θ f + θ m than θ m , which is a consequence of the estimation of θ 1 being based on D e , which is influenced by θ m and θ f .
Having established that the dispersion of a scalar in a jet is a self-similar process in § 4, we seek a similarity solution to (6.1). We express (6.1) as with C m subject to the conditions where F 0 is the source flux and λ * ≡ θ 0 + 6θ 1 is a constant that characterises the propagation speed of the front. Based on the observations reported above, λ * ≈ 1.10. Introducing the similarity variable Transformation of (6.2) according to (6.5) and (6.4) results in where Pe ≡ λ * M 1/2 m /(2αD e ) ≈ 10.99, based on our observations. The transformed boundary condition and integral constraint become Equation (6.6) has the general solution Figure 9 demonstrates that the similarity solution (6.8) agrees reasonably well with the observed data. In particular, the model reproduces the shape and the longitudinal extent of the scalar's distribution. The model solution appears to slightly overpredict the peak concentration, and the concentration in the tail of the distribution at relatively large values of λ. We attribute the latter to higher-order contributions in the dimensionless scalar flux θ m , for which (5.2) does not account (cf. the higher-order terms for dispersion in a pipe flow considered by Chatwin 1970). However, given that the model is based on a leading-order representation of the front involving two parameters, and that the specification of θ 0 was made a priori, the agreement is satisfactory.

Planar jets
Although the dispersion model described in this paper has been developed for axisymmetric jets, it is built from general principles that are equally applicable in the analysis of planar jets or plumes. Central to the model are the geometric parameters θ 0 = f 0 g 0 and θ 1 = f 0 g 1 , which correspond to the dimensionless scalar flux of the steady state and of the perturbed scalar profile respectively. Both θ 0 and θ 1 are functions of the dimensionless longitudinal velocity f 0 (η), the scalar concentration g 0 (η) and the perturbed concentration g 1 (η). Due to the fact that g 1 (η) can be determined a priori by solving an ordinary differential equation involving g 0 , the dimensionless perturbation flux θ 1 can, in principle, be determined without consulting unsteady data.
In contrast, Landel et al. (2012), while providing several analytical solutions and insights from experiments, relied on advection and dispersion parameters (K a and K d respectively), whose relationship with the underlying turbulent scalar flux w c or mean flow scalar flux w c was not made explicit. However, the validity of both the model presented by Landel et al. (2012) for dispersion in quasi-two-dimensional jets (1.2), and the framework we propose is strengthened by the fact that when our dispersion Shear-flow dispersion in turbulent jets 45 closure is applied to planar jets we find (see appendix B) which is in one-to-one correspondence with (1.2) when It is important to note, however, the fundamental physical difference between a statistically two-dimensional turbulent planar jet and the quasi-two-dimensional jet studied by Landel et al. (2012) by confining the jet to the narrow gap between two parallel plates. Although our treatment in appendix B applies equally well to each of these cases, one would expect dimensionless parameters such as the entrainment coefficient α and θ 1 to depend on the particular case considered. Not straightforward to anticipate from (1.2) is that the advection parameter K a depends on both θ 0 and θ 1 . Using (7.2), it is useful to compute the values of θ 0 and θ 1 that correspond to K a = 1.65 and K d = 0.09, as reported by Landel et al. (2012): where the values of θ 0 and θ 1 that we report in the present paper are shown in parentheses, and we have taken α = √ 2 × 0.068, based on the measurements of Landel et al. (2012) (the factor √ 2 arising due to the characteristic scales we use to define α). Whereas K a and K d were obtained by fitting to experimental data, θ 0 can be obtained theoretically as the dimensionless scalar flux arising from Gaussian profiles of w and c of equal spread. On the other hand, θ 1 appears to take different values in axisymmetric jets compared with quasi-two-dimensional jets, which is perhaps not surprising when one considers that θ 1 depends on the lateral mixing provided by turbulence and the longitudinal scaling of the velocity field. In planar jets, the dispersion closure used in (7.1) has the form (see appendix B) While it may appear inconvenient that (7.5) contains prefactors that differ from those in the equivalent expression for axisymmetric jets (5.8), the advantage is that the physical meaning of θ 0 and θ 1 is unchanged. In principle, one is therefore able to make predictions a priori about dispersion in planar/quasi-two-dimensional jets in addition to axisymmetric jets.

Conclusions
We have analysed the transport of a passive scalar in a statistically axisymmetric turbulent jet and found that the longitudinal mixing is primarily the result of shear-flow dispersion rather than turbulent transport (see, e.g. figure 8). Consequently, estimations of the amount of longitudinal mixing based exclusively on turbulence that do not account for the radial dependence of the mean flow are likely to significantly underestimate the total longitudinal mixing.
While shear-flow dispersion in turbulent jets is ostensibly different from shear-flow dispersion in pipe flows, for which the classical theory was first developed (Taylor 1953), there are several similarities. Different from pipe flow is the fact that in a jet the radial steady-state scalar concentration is non-uniform, and will, according to the lateral transport properties of the scalar, be concentrated in regions of relatively high longitudinal velocity. From an integral perspective, this property of the scalar field is characterised by a constant dimensionless flux parameter θ 0 . Similar to pipe flow is the fact that longitudinal changes in an otherwise steady-state scalar field result in a local perturbation of the scalar's radial dependence. In jets, the effect of a perturbation of amplitude L 1 on the local scalar flux can be accounted for with the parameter θ 1 . Together, θ 0 and L 1 θ 1 provide a dispersion closure in an integral model that shows good agreement with DNS data (figure 9). When applied to planar jets the approach yields a differential equation that is in one-to-one correspondence with that proposed by Landel et al. (2012) and provides additional physical insight.
In principle, shear-flow dispersion finds application in any problem in which there are lateral gradients in a velocity field, and is therefore relevant to a wide range of situations in industrial and environmental fluid mechanics. In particular, one should expect to find shear-flow dispersion in both laminar and turbulent free-shear flows, including turbulent plumes and wakes. The axisymmetric turbulent jet that we have analysed is relatively simple because it is statistically two-dimensional and self-similar at relatively large distances from the source. Consequently, the behaviour of a propagating step change in the integral concentration could be characterised completely in terms of a similarity variable λ ∝ z 2 /t. However, one expects shear-flow dispersion to also play a role in more complicated three-dimensional problems such as plumes that are influenced by a cross-wind and multiple coalescing jets and plumes. The present study indicates that unless shear-flow dispersion, arising from advection by the mean flow, is correctly understood and parameterised, longitudinal mixing in such flows will not be predicted correctly.
Outstanding among issues deserving further attention is a formulation for the precise way in which the scalar profile is perturbed as a result of a step change in the longitudinal concentration profile. Such a formulation would allow one to estimate a dispersion coefficient (cf. θ 1 in the present study) for a variety of different free-shear flows such as planar jets and plumes. A related issue for consideration is the role of higher-order perturbations to the scalar profile such as those considered in the context of pipe flow by Taylor (1954a), Gill (1967) and Chatwin (1970).
with boundary layer theory we neglect second derivatives in the longitudinal direction. Unlike jets, in pipes the steady-state concentration c m0 is independent of z, and we can replace c with C in the scalar transport equation: where κ T is the (constant) eddy diffusivity and w is the steady-state velocity field. Taylor (1953) applied a Galilean transform, defining a coordinate system that moves according to the mean longitudinal velocity or, equivalently, the velocity associated with the scalar flux: where θ 0 is constant and equal to unity for pipe flow. We assume that the radial dependence of the velocity field does not vary in the longitudinal direction, such that w = w m f 0 (η), where w m = w is the cross-sectional-averaged longitudinal velocity, η ≡ r/r m is a normalised radius and f 0 = 1. Equation (A 1) transforms according to The radial dependence of C can be decomposed into g 0 (η), which corresponds to the steady-state radial dependence of the concentration, and g 1 (η), which accounts for the leading-order perturbation from g 0 (η): C ∼ C m (g 0 + L 1 g 1 ).
(A 4) Here, the variable C m determines the amplitude of the concentration profile. For pipe flow, since the steady-state concentration is uniform, it can be assumed that g 0 = 1, without loss of generality, and we define C m such that the cross-sectional average g 1 = 0. In (A 4), L 1 determines the amplitude of the leading-order perturbation. Applying the averaging operation to (A 3) gives an equation for the mean concentration: ∂C m ∂τ + w m f 0 g 1 ∂ ∂ξ (C m L 1 ) = 0. (A 5) Subtracting (A 5) from (A 3) yields (A 6) Taylor's key assumption was that for large time (t r 2 m /κ T ) there exists an asymptotic balance between the second and fourth terms in (A 6). Noting that ∂ ξ = ∂ z , one therefore assumes that which represents a balance between longitudinal advection and turbulent transport in the radial direction. Equation (A 7) has a separable form and implies that to within a constant multiplicative factor In principle, (A 7) can be inverted to determine the form of the perturbation g 1 , and therefore the flux F m ≡ w m c m r 2 m ( f 0 g 0 + L 1 f 0 g 1 ).
(A 9) Defining θ m ≡ F m /w m c m r 2 m , θ 0 ≡ f 0 g 0 , θ 1 ≡ f 0 g 1 , (A 10a−c) (A 9) can be expressed as θ m = θ 0 + L 1 θ 1 . The parameter θ 0 is the basic steady-state dimensionless scalar flux arising from self-similar profiles of c, and θ 1 is the dimensionless scalar flux arising from a local departure from self-similarity in c.