# The turbulence boundary of a temporal jet

###### Abstract

We examine the structure of the turbulence boundary of a temporal plane jet at using statistics conditioned on the enstrophy. The data is obtained by direct numerical simulation and threshold values span 24 orders of magnitude, ranging from essentially irrotational fluid outside the jet to fully turbulent fluid in the jet core. We use two independent estimators for the local entrainment velocity based on the enstrophy budget. The data show clear evidence for the existence of a viscous superlayer (VSL) that envelopes the turbulence. The VSL is a nearly one-dimensional layer with low surface curvature. We find that both its area and viscous transport velocity adjust to the imposed rate of entrainment so that the integral entrainment flux is independent of threshold, although low-Reynolds-number effects play a role for the case under consideration. This threshold independence is consistent with the inviscid nature of the integral rate of entrainment. A theoretical model of the VSL is developed that is in reasonably good agreement with the data and predicts that the contribution of viscous transport and dissipation to interface propagation have magnitude and , respectively. We further identify a turbulent core region (TC) and a buffer region (BR) connecting the VSL and the TC. The BR grows in time and inviscid enstrophy production is important in this region. The BR shows many similarities with the turbulent-nonturbulent interface (TNTI), although the TNTI seems to extend into the TC. The average distance between the TC and the VSL, i.e. the BR thickness is about 10 Kolmogorov length scales or half a Taylor length scale, indicating that intense turbulent flow regions and viscosity-dominated regions are in close proximity.

M. van Reeuwijk and M. Holzner]MAARTENvanREEUWIJK^{†}^{†}thanks: Email address for correspondence: m.vanreeuwijk@imperial.ac.uk

and

MARKUS HOLZNER

## 1 Introduction

Turbulent entrainment is the incorporation of ambient fluid at the boundary of turbulent flows such as free shear flows or at the free stream edge of turbulent boundary layers. It is an important process in a variety of engineering and geophysical flows controlling the turbulent transfer of mass, heat and momentum (Da Silva et al. 2013; Pope 2000; Stull 1998; Thorpe 2005). A relevant, yet unresolved issue that has received renewed interest in recent years is the connection between processes that are associated with the large scale organization of the flow and processes that occur at the scale of the smallest eddies (e.g. Westerweel et al. 2005; Da Silva & Taveira 2010; Hunt et al. 2011; Philip & Marusic 2012; Wolf et al. 2012).

The integral rate at which ambient fluid is incorporated into the turbulent flow, in the following referred to as global entrainment, is independent of the small scale details of the flow, i.e. it does not depend on the viscosity or the energy dissipation mechanism. The common entrainment assumption is that the global entrainment velocity is proportional to the typical velocity inside the turbulent zone (Morton et al. 1956; Turner 1986), usually the centreline velocity. The entrainment coefficient is typically , but the value is far from universal and depends on the choice of the typical length scale , the assumed shape of the velocity profile and can also depend on the initial conditions (e.g. Redford et al. 2012).

Conversely, Corrsin emphasised a microscale perspective, in the following referred to as local entrainment, and suggested that the turbulence boundary is demarcated by a very thin viscosity-dominated laminar superlayer, whose local propagation velocity towards the non-turbulent region is determined by two parameters: the kinematic viscosity and the rate of dissipation of kinetic energy (Corrsin & Kistler 1955). Consequently, where is the Kolmogorov velocity scale. The ratio of local entrainment velocity to global entrainment velocity is given by

(1.0) |

and the Reynolds number dependence begs the question in which way the two views - local and global - are consistent. The dependence on Re seems to imply that both surface area and viscous diffusion adjust to the imposed global entrainment rate such that the small scale details of how the vorticity is transferred are somehow forgotten across interactions of eddies with a large hierarchy of sizes (Townsend 1976). By denoting the integral entrainment flux as , the global perspective suggests that , where is the surface area based on the average distance of the turbulence interface to the core of the turbulent zone. From the local perspective where is the total surface area of the contorted turbulence boundary. Equating the two expressions for results in and therefore

(1.0) |

This means that must be large to compensate a slow viscous transfer of vorticity and to cancel out the viscosity dependence (e.g. Tritton 1988; Sreenivasan et al. 1989).

Probably the simplest setting where turbulence propagates into non-turbulent fluid is the case without any mean flow, which can be realized via oscillating grid experiments, e.g. Holzner et al. (2007, 2008); Holzner & Luethi (2011). The results obtained in such a flow showed evidence for the presence of a laminar superlayer at the boundary of turbulent flow regions. In particular, the analysis supported that is indeed given by a strongly convoluted surface and accounts for a large entrainment flux with a small characteristic velocity comparable to the Kolmogorov velocity (Holzner & Luethi 2011). A similar picture, i.e. , recently emerged from the experiments in a round jet of Wolf et al. (2012). In all the experiments and simulations, the probability density functions (PDFs) of the entrainment velocity indicated that there is a large variation in entrainment velocities. In this context the term laminar superlayer is unfortunate as it suggests that the flow is layered without notable fluctuations. Therefore, this layer will be termed the viscous superlayer (VSL) in the remainder of the paper.

A somewhat different view emerges from direct numerical simulations of plane jets (Da Silva & Pereira 2008; Da Silva & Taveira 2010; Taveira & Da Silva 2013), a plane wake (Bisset et al. 2002) and experiments in a round jet (Westerweel et al. 2005, 2009), which focused on properties of the turbulent/nonturbulent interface (TNTI). The TNTI seems to be thicker than the VSL predicted by Corrsin; that is, the thickness of the TNTI is comparable to the Taylor length scale , rather than than the Kolmogorov length scale . The difference in character between the layers is also evident from the dominant physical processes: in the TNTI turbulence propagates mostly via transmission of turbulent (i.e. Reynolds type) shear stresses (Westerweel et al. 2005, 2009), whereas it is the action of viscous shear forces in Corrsin’s theory. Bisset et al. (2002) make an explicit distinction between the two layers by stating that the TNTI is a thin turbulent layer connecting non-turbulent (irrotational) and the turbulent regions of the flow, and they conjecture that the VSL forms the outer boundary of the TNTI. The aims of this paper are firstly to determine whether a VSL can be observed for a generic shear flow and secondly to study in detail the structure of the turbulence boundary.

One important factor that may partly explain the observed differences in layer properties is the method by which the interface between the nonturbulent and turbulent fluid is identified. Indeed, the interface is usually obtained by applying a threshold to a scalar field such as enstrophy (Bisset et al. 2002; Mathew & Basu 2002; Holzner et al. 2007; Da Silva & Pereira 2008) or a high Schmidt number dye (Westerweel et al. 2005, 2009). By construction, this interface is artificial because the transition between turbulent and nonturbulent fluid must occur smoothly over a finite region. The threshold value is generally chosen in a range where results are insensitive to the precise threshold value, e.g. for the conditional statistics (Bisset et al. 2002; Holzner et al. 2007; Westerweel et al. 2005, 2009). A complication in this matter is that in experiments and even in numerical data sets it is often difficult to vary thresholds over a span of several decades because of experimental (e.g. Westerweel et al. 2005; Holzner et al. 2007) or numerical (e.g. Bisset et al. 2002; Mathew & Basu 2002) noise.

In this paper we perform a systematic study of the effect on the threshold value on the entrainment velocity and related statistics. In doing so we span the entire range from essentially irrotational fluid near the turbulence boundary to fully turbulent fluid near the jet centre, which will enable us to address whether a VSL exists at the outer fringes of turbulence and how it may be related to a TNTI. Our analysis supports the existence of a VSL over a large range of thresholds ( 20 decades), a turbulent core and a smooth transition zone connecting the two that will be identified as the buffer region. The buffer region shows several features characteristic of the TNTI. The study provides insight into how the local and global turbulent entrainment are connected. We find that the VSL is a nearly one-dimensional layer with low surface curvature and both its area and viscous transport velocity adjust to the imposed rate of entrainment so that the entrainment flux is independent of threshold. We perform simulations of a temporal plane jet at which is a generic well-documented flow (Da Silva & Pereira 2008) that is attractive from the viewpoint that it has two homogeneous directions. The present paper starts out with a theoretical framework (§2) that describes the properties of the local entrainment velocity from a local and a integral perspective. A simple conceptual model is developed that is tested against the simulation data. After a brief explanation of the simulation setup (§3), the entrainment characteristics, layer structure and geometry of the turbulence boundary are presented (§4). In §5, the relation between the identified layer structure and the TNTI is explored, as well as the mechanism by which entrainment takes place. Concluding remarks are made in §6. In a companion paper the influence of mean shear and Reynolds number are analysed.

## 2 Theory

This section revisits the various definitions of entrainment velocity, the determination of the local entrainment velocity based on enstrophy budgets from a local perspective and complements it with an integral approach to the problem. Thereafter a simplified model is set out with predictions for the enstrophy transport across the VSL.

### 2.1 Definitions of entrainment velocity

One of the subtleties of turbulent entrainment is that there are several definitions for the entrainment velocity in use (Turner 1986; Hunt et al. 1983). The most common definition is the rate at which fluid flows into the turbulent zone across its boundary, commonly denoted . For the temporal jet, the specific flux is constant and therefore (see also §3). A second definition is the rate at which the edge of a turbulent flow spreads outwards, i.e. the boundary entrainment rate . This is the quantity that is denoted in this work by . A third definition is which can be interpreted as the speed of the interface relative to the mean fluid velocity. Note that and are global (macroscale) quantities measured in a fixed coordinate system (sometimes called laboratory coordinates, e.g. Westerweel et al. 2009), whereas uses a coordinate system relative to the mean flow. In the next section we will also define local (microscale) entrainment velocities that, similar to , represent the local interface propagation velocity relative to the local fluid velocity.

### 2.2 Local entrainment velocity: classical approach

We differentiate between turbulent and non-turbulent flow regions by using a threshold on the enstrophy , where is a component of vorticity. This defines a bounding surface separating the two regions. In a Lagrangian frame moving with the iso-enstrophy surface the isolevel will by definition remain constant and this property can be used to derive an expression for the propagation velocity (Holzner & Luethi 2011). We write the velocity of an isosurface element, , as a sum of fluid velocity, , and velocity of the area element relative to the fluid, , that is, . The total change of in a frame of reference moving with an enstrophy isosurface element is then given by

(2.0) |

where the lower case is the total derivative following a surface element, and the upper case is the material derivative which follows a fluid element. By defining a surface normal as and the normal relative velocity component , we obtain

(2.0) |

Substituting the enstrophy balance equation

(2.0) |

into (2.2) and averaging over the isosurface , we obtain an expression for the average entrainment velocity :

(2.0) |

where

### 2.3 Local entrainment velocity: integral approach

An alternative to the local approach described in 2.2 is to integrate the enstrophy equation (2.2) over a time-dependent domain which has boundary velocity and use the Reynolds transport theorem, resulting in

As the surface normal points into the turbulent region, the appropriate volume under consideration comprises the irrotational region. We can formalise this by defining a control volume where is the Heaviside function and is an enstrophy threshold. As is then by definition constant on the surface , the equation above simplifies to

(2.0) |

Using the Gauss divergence theorem and making use of the fact that for an incompressible fluid, the entrainment flux can be expressed as

(2.0) |

Introducing an average over the volume , (2.0) can be rewritten as

Now, because spans the entire nonturbulent region, it is expected that , and therefore

(2.0) |

where

### 2.4 A model for enstrophy transport in the viscous superlayer

In the VSL, the evolution of enstrophy is governed by molecular processes (Corrsin & Kistler 1955; Holzner & Luethi 2011), i.e. . Assuming that the local curvature is small on average and multiplying by , Eq. (2.2) then becomes

(2.0) |

Integrating this expression and using that at , both and , we obtain

(2.0) |

Assuming that is constant in the VSL, the square of the solution to (2.0) is

(2.0) |

where is a reference value for enstrophy. Hence, the enstrophy is expected to drop off exponentially in the VSL, provided that is constant (see §4.2).

The model solution (2.0) can be used to predict the magnitude of the enstrophy transport terms. For the local approach we expect

Entirely consistently, we expect for the integral approach that

## 3 Simulations

The start situation for a temporal plane jet is a fluid layer that is quiescent except for a thin region where the streamwise velocity is nonzero, and is homogeneous in the two other directions and . Here, is the initial jet width. It follows from continuity that the volume flux remains constant for this flow throughout the jet’s transition to turbulence and subsequent growth due to turbulent entrainment. Assuming that the Reynolds number , the only relevant parameters are the initial volume flux and time which suggests self-similar behaviour, with the jet halfwidth and centreline velocity scaling as and , respectively.

The simulation domain is a cuboid of size , which is three times larger in all directions than the domain used in Da Silva & Pereira (2008). The larger domain facilitates much longer simulations, thereby allowing not only the first moments but also the turbulence to reach an equilibrium and has the added advantage of improved statistics because of the larger area spanned by the two homogeneous directions. Periodic boundary conditions are employed in the two homogeneous directions and , and free-slip boundary conditions are imposed at .

The simulation considered here is for . The resolution of the simulation is which is sufficient to ensure that all active length scales in the turbulence are fully resolved. We define a reference time-scale and simulate for . All statistics before are discarded to ensure that the turbulence has time to reach a dynamic equilibrium.

Following Da Silva & Métais (2002); Da Silva & Pereira (2008), we use the initial condition

(3.0) |

where is chosen such that . We set (Da Silva & Pereira 2008) and seed the initial condition with uniform random noise with an enstrophy level that is 8 percent of the maximum average enstrophy [which is ]. Note that the perturbation amplitude in terms of the velocity is only about one percent. This facilitates a rapid transition to turbulence and we note that the enstrophy levels after the transition to turbulence far exceed the noise levels. The code used for direct numerical simulation is based on fully conservative second order finite difference operators in space (Verstappen & Veldman 2003) and uses an adaptive second-order Adams-Bashforth time integration scheme. The advantage of the spatial discretisation used is that the numerics are free of numerical diffusion whilst still satisfying volume and momentum conservation. More details can be found in van Reeuwijk et al. (2008).

As the statistics shown in the next section depend heavily on budgets for the enstrophy, special care is taken to ensure that the budgets are calculated consistently with the numerics used. To achieve this, a mimetic (Hyman & Shashkov 1997; van Reeuwijk 2011) curl operator is defined such that it satisfies the identity up to machine precision, where denotes pressure. In order to ensure calculation fully compatible with the numerical method, we do not manually discretise (2.2), but instead make use of the following identities

(3.0) | |||

(3.0) |

which are then also enforced on the discrete level. Taking (3.0) as an example, one can calculate the first term on the right hand side directly by using the routine for scalar diffusion on ; the second term can be calculated by taking the discrete curl of the viscous term in the momentum equation, and then taking the scalar product of the result with the vorticity components.

As the temporal jet has a nonzero mean velocity in the -direction, it is important to ensure that the identity is also satisfied on the discrete level. Indeed, this identity was used explicitly to derive (2.0). This can only be achieved if the thresholding algorithm identifies entire cells to be either inside or outside the turbulent region. Indeed, we found that if we used trilinear interpolation to construct an isosurface - which is in principle a better representation - the various interpolations required could lead to very significant deviations in the calculated entrainment velocity.

## 4 Results

### 4.1 Bulk flow properties

The time evolution of the enstrophy levels in the jet are shown in Fig. 1. Here, is a reference enstrophy level. Fig. 1(top left) shows the initial condition and the low-amplitude noise. At , the flow has fully transitioned to turbulence and high enstrophy levels can be observed in the jet that very rapidly drop off near the jet edge. As time progresses further, the enstrophy levels decrease and spatial scales can be seen to grow.

Fig. 2(a) demonstrates that the grid resolution is appropriate for the problem under consideration. Shown is the grid spacing normalised by the Kolmogorov length scale based on the centreline dissipation rate . The overbar denotes averaging over the two homogeneous directions and over 10 . The dissipation rate has its maximum at and dissipation rates will be much lower at the jet boundary, which implies that the simulation is even better resolved there (dashed and dash-dotted lines). As can be seen, the simulation becomes better resolved in time because , which can be inferred by using , as is confirmed in Fig. 4(f).

Fig. 2(b) shows the evolution of the Taylor Reynolds number, defined as , where is the centreline turbulent kinetic energy and is the Taylor microscale (e.g. Tennekes & Lumley 1972, pp. 67-68). As judged from , the turbulence reaches an equilibrium value of about after .

Energy density spectra for the plane , averaged over shells of wave number and a time interval of are shown in Figure 3. Fifteen spectra are shown for and the collapse demonstrates the self-similarity of the flow under consideration; even though the spectra change in time, the normalisation with and cause a full collapse of the data. The spectra indicate that there is a range of active scales and that there is a small separation of scales as is evident from the formation of a spectrum (left) and a comparison of the peak in the energy density (left) and dissipation spectrum (right). Note that the dissipation spectrum peaks at , once more indicating that the simulation is fully resolved.

As the velocity profile is symmetric around , the jet half-width was inferred from the average of the values of for which and . For all profiles shown, use has been made of the symmetry (or anti-symmetry) in the profile to further improve the statistical accuracy. Shown in Fig. 4(a) is the scaling of with time, which as expected from dimensional arguments is linear; the red dashed line is a linear fit. The normalised mean velocity and momentum flux are shown in Fig 4(b,c). These profiles were scaled and then further averaged over four contiguous time-intervals (an effective average over 40 ). The profiles are convincingly self-similar, although the profile of shows more variability than because it is a second order moment. Also shown in Fig 4(b) is data from laboratory experiments of a plane jet by Gutmark & Wygnansky (1976) (red upward triangles) and Ramparian & Chandrasekhara (1985) (blue downward triangles) as well as the numerical simulations of a temporal plane jet by Da Silva & Pereira (2008) (green circles). Good agreement can be observed. Self-similarity of turbulent quantities is demonstrated in Figs 4(d-f). As discussed earlier, the balance between turbulence production and dissipation suggests that (Fig. 4(f)), which in turn implies that (Fig. 4(d)). The profiles for the turbulence kinetic energy (TKE) and the dissipation rate (Figs 4(e, f)) show once more a reasonably good collapse, with peaking at (where shear-production is maximal) and the dissipation rate peaking at the centreline. We observed variations in the growth rate of between different simulations upon variation of the initial conditions, despite a convincing self-similar behaviour in all of them. This may point to non-universal self-similar behaviour. Indeed, Redford et al. (2012) showed through simulations of an axisymmetric temporal wake that differences in the initial conditions can influence growth rates (and therefore entrainment coefficients) for extended periods. During such periods, the flow developed in a self-similar fashion but was not universal; only on very large timescales was universal behaviour observed. The non-universal self-similarity may explain the slightly higher turbulence levels in the current simulation compared to those observed by others (Fig. 4(e)).

### 4.2 Entrainment velocity and budgets

A cross-section of the enstrophy field at constant is shown in Fig. 5, together with enstrophy isolines at and . At , the turbulent/nonturbulent interface is highly contorted and has ’holes’, whereas the lower thresholds do not have such holes. Note that what may appear as holes on the figure (i.e. in 2D) for low tresholds are fluid portions that are connected in 3D to the outer irrotational region, while at higher thresholds one also finds islands of low vorticity disconnected from the outer region. What is striking is that the enstrophy levels at and are in close proximity to the high enstrophy regions, indicating a very quick drop-off of enstrophy levels near the turbulence boundary. At these low threshold values, the surface remains contorted because of the large scale vortices distorting the flow but the enstrophy isosurfaces appear to form nearly one-dimensional layers with relatively small curvature.

The instantaneous budgets of enstrophy were calculated for the entire 3-D field every 5 and then used to calculate the terms in (2.2), (2.0) for 37 thresholds in the range . Simultaneously, the volume was recorded for each of the threshold values, which enabled an independent calculation of the propagation velocity using (2.0). Shown is the data for .

First, we show that the calculated local entrainment velocities and correspond to the actual local entrainment velocity . To this end, the terms comprising the local entrainment velocity (2.2) and (2.0) normalised by the directly measured volume-based entrainment velocity are presented in Fig. 6. For , the predicted propagation velocities (squares, Fig. 6(a)) and (squares, Fig. 6(b)) match the actual propagation velocity excellently. A small systematic error can be discerned in the calculation of , as the calculated propagation velocity shows a small but systematic trend in . This systematic trend is not observed in , although the prediction is slighty lower than . The poor predictions for are not associated with and but are due to an insufficient temporal sampling frequency creating large errors in the calculation of (Eq. 2.0); this occurs because the high enstrophy regions tend to shrink and expand rapidly on a timescale shorter than .

Fig. 6(a) shows convincing evidence of the existence of a viscous superlayer (VSL). Indeed, in the VSL the inviscid contribution (circles) does not play a role (Holzner & Luethi 2011) and we observe that this is the case for . For , the inviscid terms increase very rapidly. Also shown in Fig. 6(a) are the theoretical predictions from §2.4, namely and (both displayed by dotted lines). The predictions are in good agreement, although the observed magnitudes of (downward triangles) and (upward triangles) are a bit larger than predicted, and a dependence on is discernable.

The budget of is similar to that of : the temporal and inviscid contributions, (diamonds) and (circles) respectively, are negligible and the propagation of the enstrophy isosurface in the VSL is caused by viscous effects. For this indicator, the contributions and seem to become fully independent of the threshold level below .

One of the main assumptions made in the derivation of the theoretical model was that the curvature of the isosurface is small. The reasonably good agreement with the theoretical model in Fig. 5 supports this assumption but by using Eq. 2.0 it can be validated explicitly. Figs 6(c,d) show the contribution to by curvature (red crosses) and diffusive transport in the direction of the surface normal (blue triangles). For both as , the curvature term becomes negligible for . This provides confirmation that the theoretical model in §2.4 is a reasonable description of the processes governing the VSL.

Having established that both estimates of the local entrainment velocity and are in good agreement with the actual entrainment velocity , we study the dependence of the entrainment velocity on the threshold value . In Fig. 7(a) we show the dependence of , and on , normalised by the Kolmogorov velocity scale . The vertical error bars denote the variation over the entire time interval . There is a clear dependence of on the enstrophy threshold: isosurfaces for the very low enstrophy thresholds propagate faster than those with higher thresholds. Indeed, is nearly twice as high at than at . Hence, although is of the same order of magnitude as , the dependence on suggests that it is not merely the viscosity and dissipation rate that determine the propagation velocity in the VSL. This may be a low-Reynolds number effect: as the Reynolds number increases, the VSL will become thinner relative to the integral scale and therefore the enstrophy levels in the VSL will drop off quicker, cf. (2.0), leaving less opportunity for variation in . It should also be noted that the surface is not smooth but follows the grid (cf. §3) to ensure conservation properties. Further work is required to settle this issue.

Another striking feature is that becomes zero around and positive for . Hence, high enstrophy regions move inwards towards the jet centre, low enstrophy regions move outwards and there exists an isosurface that separates the shrinking and expanding regions. The movement of high enstrophy regions towards the jet centre can be explained by using the relation between enstrophy and the dissipation rate which is valid for isotropic and homogeneous turbulence . Using the self-similarity of it follows that . Hence, if one would assume that has a self-similar profile and pick a reference threshold that one follows in time, it would be seen to move inwards towards the jet core. This applies to enstrophy levels where turbulence is developed and is positive. Towards the VSL, at low levels, viscous transport is diffusing enstrophy outwards and is negative. Note that the gradient of over in Fig. 7(a) is always negative meaning that the enstrophy profile is flattening over time.

As macroscale entrainment is independent of molecular processes, it is expected that the global entrainment velocity is independent of threshold. As mentioned in §2, for the case under consideration corresponds to the boundary velocity in the terminology of Turner (1986) and is defined as where . The relation to (Eq. 2.0) is therefore . Fig. 7(b) shows the global entrainment velocity normalised by the centreline velocity amplitude . The standard entrainment assumption (e.g. Turner 1986) is and the ratio plotted is the entrainment coefficient . As can be seen, is practically constant as a function of although there is a small systematic trend. When is independent of , it implies that the entrainment flux is constant in the VSL; indeed, by using it immediately follows that is constant if is independent of . This result indicates that to first order , and only if is independent of do we expect that independently of , i.e. the classical view advocated by Corrsin.

The small trend discernable in Fig. 7(b) is a low-Reynolds effect associated with the position of the interface. Indeed, selfsimilarity implies that , where is the universal velocity profile and the selfsimilarity variable. The entrainment assumption is , where is defined as and in our case is the half-velocity width. As mentioned in §1, other definitions of will result in different values for . This is straightforward to see by fixing at a value differing from unity. Indeed, for an alternative width , we obtain

(4.0) |

indicating that the effective entrainment rate is . Hence, if there is a dependence of the average interface position on this will create a trend in . In Fig. 7(c), we have plotted , where and is the mean -position of the isosurface. As can be seen, the value of is constant for both the directly measured entraiment velocity and the calculated entrainment velocity . There is a small downward trend for the value of as calculated from which is due to the small systematic error discussed earlier. The dependence of on will no longer play a role at very high Reynolds numbers because the VSL will become so thin that will become independent of .

In summary, the dependence of the variation of on can be explained by two independent mechanisms: 1) the dependence of the surface area on ; and 2) the finite thickness of the turbulence boundary as scaled on the jet thickness . These are likely to be low-Reynolds-number effects and it is expected that in the limit of , the classical Corrsin viewpoint will be recovered. Further work is required to verify this hypothesis.

### 4.3 The structure of the turbulence boundary

Having access to the entrainment velocity as a function of the threshold value allows one to explore the structure of the turbulence boundary. There are two distinguishing enstrophy threshold values to consider at any time, namely (i) the enstrophy level at which and (ii) the enstrophy level for which enstrophy production becomes negligible, diagnosed by the criterion , where is a small value. As discussed in the previous section, the level is a threshold that separates expanding regions ( from shrinking regions () and the threshold demarcates the start of the VSL. This suggests a layer structure as shown in Fig. 8. Away from the jet boundary, the flow is non-turbulent and irrotational. Moving closer to the turbulence boundary, one enters the VSL. The transition location is arbitrary and would depend on a choice of threshold. The VSL extends up to the location where inviscid terms start playing a role. Note that, similar to the viscous sublayer in a wall-bounded flow, the VSL can be classified neither as turbulent nor as laminar because viscous effects are dominant whilst there are significant fluctuations in the layer due to external influences. We define the turbulent core (TC) to be the region for which , and define a buffer region (BR) to be the zone between the VSL and the TC, which is still expanding but for which inviscid terms are important. This term was chosen to emphasise the resemblance with the buffer layer of a wall-bounded turbulent flow that also couples two regions in which different processes dominate (viscous and inertial effects in the case of a wall-bounded flow).

The evolution of the layer structure is plotted as a function of time in Fig. 9 for (dashed lines) and (solid lines). A value was adopted to demarcate the onset of the VSL. Although the exact values obtained by the two estimates differ, the trends are consistent. The three regions are clearly visible in Fig. 9. As time progresses, the threshold level for which moves to lower and lower values, in accordance with the decay expected from self-similarity. The threshold level demarcating the beginning of the VSL remains approximately constant at for , although significant fluctuations can be observed for and the threshold values differ by a factor in the range where the VSL resides. Hence, it is impossible to infer unambiguously whether the onset of the VSL occurs at a fixed threshold value or not.

### 4.4 Geometry: relating to a distance

The dependence of volume and surface area on the threshold value can be used to obtain information about the average distance from one enstrophy isosurface to the next, thereby getting an impression of the distance between different regions of the flow. Making use of the observation that curvature is low at low thresholds, we can define a average distance , which is related to the volume for which and the surface area as

(4.0) |

Note that the sum of and is exactly half of the simulation domain volume. Introducing and as the difference in volume and distance respectively between two adjacent enstrophy thresholds and , the average distance can be approximated by . We set at the start of the VSL at .

The dependence of on is shown for in Fig. 10(a). The collapse of the profiles for different times shows that the result is quite robust. It is striking how close essentially irrotational regions are located to regions with very high enstrophy levels on average. Indeed, the BR, which connects the TC and the VSL, is on average about or thick. Moreover, the full 24 decades in enstrophy levels are separated by or . Also shown is the analytical solution (2.0) (dotted line), expressed as with . It is evident that (2.0) is a reasonable approximation for , which explains the agreement of the simulation data with the predictions for and . However, as can be seen, the shape of the profile deviates from a straight line in a semilog plot because of the influence of surface area on . Indeed, the model assumes that is constant and this has been clearly shown not to be the case for this flow at the given Reynolds number 5000. Fig. 10(b) shows the normalized distance for a number of thresholds as a function time. After , the distances become nearly constant, indicating that the jet edge has reached a dynamic equilibrium and the distances between isosurfaces scale with .

## 5 Discussion

### 5.1 Relation with the TNTI

In the previous section we explored the properties of the turbulence boundary using a large range of enstrophy thresholds and identified a viscous superlayer, a buffer region and a turbulent core. In this section, the relation to the TNTI is established.

In their studies of the TNTI, Bisset et al. (2002); Da Silva & Pereira (2008); Mathew & Basu (2002) used an enstrophy-based threshold of . Based on Fig. 9, this would correspond to an interface located roughly on the boundary between the buffer region and the turbulent core. However, a direct comparison of the value of the thresholds might not be the best way to ascertain where the TNTI resides. Indeed, the enstrophy levels are both time-dependent and Reynolds-number dependent. This can be made explicit by using the relation valid for homogeneous turbulence, using (Fig. 4) and substituting the definition for Re, which yields

(5.0) |

Since for this particular flow , the relation above implies that and enstrophy levels thence depend both on time and on the Reynolds number. As the simulations performed in the present work are much longer than is usual, the enstrophy levels will be different than those reported by others. Indeed, if one compares the isosurface of in Fig. 5 (at ) with that shown in Da Silva & Pereira (2008) (at ) one observes that Fig. 5 has many more ”holes”. This suggests that an enstrophy threshold relative to the levels inside the turbulent core or taking the self-similarity into account is preferable to ensure a consistent interface detection over extended periods of time.

A characteristic feature of the TNTI is that when moving into the turbulent layer from the detected interface, the enstrophy quickly increases, peaks and then saturates at a fixed value slightly lower than the peak (Bisset et al. 2002; Westerweel et al. 2005; Da Silva & Pereira 2008). In Fig. 11(a), the enstrophy is plotted as a function of , for the period . In principle, this plot shows the same information as Fig. 10(a) but the axes are now linear. The viscous superlayer (VSL), buffer region (BR) and turbulent core (TC) are shown for convenience; the region for which moves outward as time progresses and is denoted by a grey area. Consistent with Bisset et al. (2002); Da Silva & Pereira (2008); Mathew & Basu (2002), the enstrophy can be seen to increase very rapidly in the transition from the BR to the TC. There is no plateau in enstrophy because the statistics presented here were obtained by conditioning on the enstrophy, not on the distance to the interface.

Another characteristic feature of the TNTI is a rapid change in the streamwise momentum (Mathew & Basu 2002; Bisset et al. 2002; Westerweel et al. 2005; Da Silva & Pereira 2008). In Fig. 11(b), the conditioned streamwise velocity is plotted as a function of . It is clear that the chosen origin is not ideal as there is no full data collapse, but it is evident that within the BR, increases rapidly. From figure 11, it can be concluded that the TNTI likely comprises part of the BR and part of the TC. The TNTI does not contain the VSL and we thus conclude that the conjecture made by Bisset et al. (2002) is correct - the VSL forms the outer boundary of the TNTI. The two are dynamically different and will consequently behave differently.

### 5.2 What makes the interface propagate?

One may speculate on the mechanism by which the turbulence boundary moves outward. Mathew & Basu (2002) present an argument on how nibbling by small-scale eddies on the Kolmogorov microscale is compatible with an inviscid macroscale entrainment process using the fractal properties of the TNTI. Conversely, Hunt et al. (2008) present an argument that larger scales of the order of the thickness of the interfacial shear layer are responsible for the net movement of the interface using a conceptual model of an eddy impinging onto a shear layer. This paper shows that the VSL is a very thin and almost one-dimensional layer governed by molecular processes that envelopes the turbulence. Below we argue that the VSL is maintained by a balance between molecular processes and the creation of a large surface area by motions over a range of scales, thereby creating a dynamic equilibrium with associated entrainment velocity of .

Indeed, the net mean motion of the VSL can be inferred from the integral scale entrainment flux , which can be expressed in turbulence quantities (as characterized by ) by because of self-similarity (Fig. 4). Due to the fractal nature of the interface, the surface area on the Kolmogorov scale is given by , where is the fractal dimension (Sreenivasan 1991). Consequently,

(5.0) |

where an equilibrium across scales of the form was assumed in the second step. The argument above suggests in the VSL , consistent with Corrsin’s argument, previous (Holzner et al. 2008; Holzner & Luethi 2011; Wolf et al. 2012) and the present work. The fractal geometry argument suggests that is convoluted over a range of scales from large to small. Referring back to Fig. 6(c,d), we show that on average the curvature of the VSL is low. This would exclude vigorous mixing on the Kolmogorov scales as a dominant process as this would create very strongly curved surfaces. Hence, it seems more plausible that motions on larger scales are more significant in close proximity to the VSL, as the vorticity is oriented and stretched tangentially to the VSL remaining correlated over larger distances, while diffusing viscously in normal direction along which it decays very sharply (Holzner & Luethi 2011). The surface area will adjust to the molecular processes governing the VSL by stretching until the product balances the inviscid entrainment volume flux .

## 6 Concluding remarks

In this paper we studied the structure of the turbulence boundary of a temporal plane jet. We find convincing evidence for the existence of the viscous superlayer (VSL). Consistent with earlier work, we find that inertial processes are negligible in the VSL. The VSL is discernible for nearly twenty orders of magnitude in enstrophy threshold. Taking into account the entire range of thresholds explored, one may attribute a thickness on the order of 15 or to the VSL, which is present for roughly . However, it should be emphasized that there is no natural threshold to define the boundary between the irrotational fluid outside the jet and the VSL. The lower this threshold is, the thicker the VSL will seem. Holzner & Luethi (2011) quantify the VSL thickness by defining it as , which in the light of the conceptual model for the VSL, Eq. (2.0), corresponds to the e-folding length which is . The simple theoretical model derived in §2.4 is in good agreement with the data and shows that the contribution of the viscous transport term amounts to and the viscous destruction term to .

The simulations support the classical assumption that global entrainment is independent of molecular processes, which was clear from the fact that the entrainment flux was practically independent of the threshold in the VSL. The local entrainment velocity was of the same order of magnitude as , although there was a dependence of the value of on the enstrophy threshold. This suggests that Corrsin’s dimensional arguments may need modification for the moderate Reynolds number under consideration here. Indeed, in the VSL both viscous transport velocity and surface area adjust to the imposed global rate so that the product is practically independent of . As the VSL becomes less contorted when moving out of the turbulent region needs to become larger to maintain a constant entrainment flux. The small dependence of on could be explained by taking into account that the average position of the interface has a non-negligible variation across the turbulence boundary.

Three regions were observed for the flow under consideration. For roughly , we observed the VSL which is characterised by a propagation velocity only depending on viscous processes. The turbulent core region (TC) was categorised as that region of the flow for which , i.e. regions that become smaller in time and that are responsible for the behaviour. The enstrophy threshold at which this region was observed was a decreasing function of time. In between the VSL and the TC there is a buffer region (BR) providing a smooth transition from the VSL and the TC. The buffer region shows many similarities with the TNTI, although the TNTI seems to extend into the TC. The core region is on average about or away from the VSL, this defines the thickness of the BR. The current work suggests that the VSL forms the outer boundary of the TNTI, confirming the conjecture made by Bisset et al. (2002). Further work, particularly a study of the Re dependence, is necessary to clarify the similarities and differences between the BR and TNTI.

Our analysis spans 24 orders of magnitude in enstrophy levels, generalizing previous approaches that are mostly based on a single threshold value. This systematic approach hence allowed us to overcome the degree of arbitrariness associated with the choice of a single threshold. Other approaches (e.g. automatic approaches) to choose an appropriate threshold for the identification of the turbulence boundary may be possible. As long as such a method is not in place, we advocate here that the analysis of a complete span of thresholds is advisable, especially because it allows separating between physically distinct regions that constitute the turbulence boundary, namely TC, BR and VSL.

## Acknowledgements

The computational resources for this work were provided by the Imperial College HPC facilities. Both authors benefited greatly from an inspiring MULTIFLOW workshop on the Turbulent/Non-turbulent Interface which took place in October 2012 in Madrid. M. H. acknowledges financial support from the Swiss National Science Foundation (SNF) under grant number PBEZP2-127831 and from the European Commission under the Marie Curie Intra-European Fellowship, project number 272886.

- Bisset et al. (2002) Bisset, D. K., Hunt, J. C. R. & Rogers, M. M. 2002 The turbulent/non-turbulent interface bounding a far wake. J. Fluid Mech. 451, 383–410.
- Corrsin & Kistler (1955) Corrsin, S. & Kistler, A.L 1955 Free stream boundaries of turbulent flows. Tech. Rep. 1244. NACA.
- Da Silva et al. (2013) Da Silva, C. B., Hunt, J. C. R., Eames, I. & Westerweel, J. 2013 Interfacial layers between regions of different turbulence intensity. Annu. Rev. Fluid Mech. 46, 457–490.
- Da Silva & Métais (2002) Da Silva, C. B. & Métais, O. 2002 On the influence of coherent structures upon interscale interactions in turbulent plane jets. J. Fluid Mech. 473, 103–145.
- Da Silva & Pereira (2008) Da Silva, C. B. & Pereira, J. C. F. 2008 Invariants of the velocity-gradient, rate-of-strain, and rate-of-rotation tensors across the turbulent/nonturbulent interface in jets. Phys. Fluids 20 (5), 055101.
- Da Silva & Taveira (2010) Da Silva, C. B. & Taveira, R. P. 2010 The thickness of the turbulent/nonturbulent interface is equal to the radius of the large vorticity structures near the edge of the shear layer. Phys. Fluids 22 (5), 121702.
- Gutmark & Wygnansky (1976) Gutmark, E. & Wygnansky, I. 1976 The planar turbulent jet. J. Fluid Mech. 73, 465.
- Holzner et al. (2007) Holzner, M., Liberzon, A., Nikitin, N., Kinzelbach, W. & Tsinober, A. 2007 Small-scale aspects of flows in proximity of the turbulent/nonturbulent interface. Phys. Fluids 19 (7), 071702.
- Holzner et al. (2008) Holzner, M., Liberzon, A., Nikitin, N., Luethi, B., Kinzelbach, W. & Tsinober, A. 2008 A Lagrangian investigation of the small-scale features of turbulent entrainment through particle tracking and direct numerical simulation. J. Fluid Mech. 598, 465–475.
- Holzner & Luethi (2011) Holzner, M. & Luethi, B. 2011 Laminar superlayer at the turbulence boundary. Phys. Rev. Lett. 106 (13), 134503.
- Hunt et al. (2011) Hunt, J. C. R., Eames, I., Da Silva, C. B. & Westerweel, J. 2011 Interfaces and inhomogeneous turbulence. Phil. Trans. R. Soc. A 369, 811–832.
- Hunt et al. (2008) Hunt, J. C. R., Eames, I. & Westerweel, J. 2008 Vortical interactions with interfacial shear layers. In IUTAM Symposium on Computational Physics and New Perspectives in Turbulence, IUTAM Bookseries, vol. 4, pp. 331–338.
- Hunt et al. (1983) Hunt, J. C. R., Rottman, J.W. & Britter, R.E. 1983 Some physical processes involved in the dispersion of dense gases. In Proc. UITAM Symp. on Atmospheric dispersion of heavy gases and small particles (ed. G. Ooms & H. Tennekes), pp. 361–395. Springer.
- Hyman & Shashkov (1997) Hyman, J.M. & Shashkov, M. 1997 Natural discretizations for the divergence, gradient, and curl on logically rectangular grids. Comput. Math Appl. 33, 81–104.
- Mathew & Basu (2002) Mathew, J. & Basu, A.J. 2002 Some characterstics of entrainment at a cylindrical turbulence boundary. Phys. Fluids 14, 2065.
- Morton et al. (1956) Morton, B. R., Taylor, G. I. & Turner, J. S. 1956 Turbulent gravitional convection from maintained and instantaneous sources. Proc. R. Soc. Lond. 234, 1–23.
- Philip & Marusic (2012) Philip, J. & Marusic, I. 2012 Large-scale eddies and their role in entrainment in turbulent jets and wakes. Phys. Fluids 24, 055108.
- Pope (2000) Pope, S. B. 2000 Turbulent flows. Cambridge University Press.
- Ramparian & Chandrasekhara (1985) Ramparian, R. & Chandrasekhara, M.S. 1985 LDA measurements in plane turbulent jets. ASME J. Fluids Eng. 107, 264.
- Redford et al. (2012) Redford, J. A., Castro, I. P. & Coleman, G. N. 2012 On the universality of turbulent axisymmetric wakes. J. Fluid Mech. 710, 419–452.
- van Reeuwijk (2011) van Reeuwijk, M. 2011 A mimetic mass, momentum and energy conserving discretization for the shallow water equations. Comp. Fluids 46, 411Â416.
- van Reeuwijk et al. (2008) van Reeuwijk, M., Jonker, H. J. J. & Hanjalić, K. 2008 Wind and boundary layers in Rayleigh-Bénard convection. i. analysis and modelling. Phys. Rev. E 77, 036311.
- Sreenivasan (1991) Sreenivasan, K. R. 1991 Fractals and multifractals in turbulence. Annu. Rev. Fluid Mech. 23, 539–600.
- Sreenivasan et al. (1989) Sreenivasan, K. R., Ramshankar, R. & Meneveau, C. 1989 Mixing, entrainment and fractal dimensions of surfaces in turbulent flows. Proceedings of the Royal Society of London Series A-mathematical Physical and Engineering Sciences 421 (1860), 79–&.
- Stull (1998) Stull 1998 An introduction to boundary layer meteorology. Kluwer Academic Publishers.
- Taveira & Da Silva (2013) Taveira, R. P. & Da Silva, C. B. 2013 Kinetic energy budgets near the turbulent/nonturbulent interface in jets. Phys. Fluids 25, 015114.
- Tennekes & Lumley (1972) Tennekes, H. & Lumley, J. L. 1972 A First course in Turbulence. MIT Press.
- Thorpe (2005) Thorpe, S. A. 2005 The turbulent ocean. Cambridge University Press.
- Townsend (1976) Townsend, A. A. 1976 The structure of turbulent shear flow. Cambridge University Press.
- Tritton (1988) Tritton, D. J. 1988 Physical fluid dynamics. Calendon, Oxford.
- Turner (1986) Turner, J. S. 1986 Turbulent entrainment: the development of the entrainment assumption, and its application to geophysical flows. J. Fluid Mech. 173, 431–471.
- Verstappen & Veldman (2003) Verstappen, R. W. C. P. & Veldman, A. E. P. 2003 Symmetry-preserving discretization of turbulent flow. J. Comput. Phys. 187 (1), 343–368.
- Westerweel et al. (2005) Westerweel, J., Fukushima, C., Pedersen, J. M. & Hunt, J. C. R. 2005 Mechanics of the turbulent-nonturbulent interface of a jet. Phys. Rev. Lett. 95 (17), 174501.
- Westerweel et al. (2009) Westerweel, J., Fukushima, C., Pedersen, J. M. & Hunt, J. C. R. 2009 Momentum and scalar transport at the turbulent/non-turbulent interface of a jet. J. Fluid Mech. 631, 199–230.
- Wolf et al. (2012) Wolf, M., Lüthi, B., Holzner, M., Krug, D., Kinzelbach, W. & Tsinober, A. 2012 Investigations on the local entrainment velocity in a turbulent jet. Phys. Fluids 24 (10), Attached publisher’s note: Erratum: ’Investigations on the local entrainment velocity in a turbulent jet’ [Phys. Fluids 24, 105110 (2012)], Phys. Fluids 25, 019901 (2013).