1 Introduction
In unbounded (e.g. jets, wakes, mixing layers) and semi-bounded (e.g. boundary layers) turbulent flows, a sharp and highly contorted interface separates the turbulent flow region from the non-turbulent ambient flow (Corrsin & Kistler Reference Corrsin and Kistler1954; Dimotakis Reference Dimotakis2000; da Silva et al. Reference da Silva, Hunt, Eames and Westerweel2014). Across this so-called turbulent/non-turbulent interface (TNTI), surrounding irrotational fluid is continuously incorporated into the turbulent flow. This process, known as turbulent entrainment, is of importance in many practical applications, in that it governs the spreading rate, mixing and reactions in a wide range of industrial and environmental flows (Simpson Reference Simpson1999; Murthy Reference Murthy2013; Davidson Reference Davidson2015).
Commonly, the TNTI is identified through a threshold on a scalar quantity, such as vorticity magnitude or enstrophy (Bisset, Hunt & Rogers Reference Bisset, Hunt and Rogers2002), turbulent kinetic energy (Holzner et al. Reference Holzner, Liberzon, Guala, Tsinober and Kinzelbach2006) or passive scalars (Westerweel et al. Reference Westerweel, Fukushima, Pedersen and Hunt2005). From a local perspective, the entrained volume flux can be expressed as the product of 
$\langle v_{n}\rangle$, the average ‘local’ entrainment velocity, where 
$\langle \cdot \rangle$ denotes an average over the surface area of the TNTI, and 
$A_{\unicode[STIX]{x1D702}}$, the surface area of the TNTI (Sreenivasan, Ramshankar & Meneveau Reference Sreenivasan, Ramshankar and Meneveau1989; Mathew & Basu Reference Mathew and Basu2002). To date, it is widely accepted that 
$A_{\unicode[STIX]{x1D702}}$ has a fractal shape (Sreenivasan et al. Reference Sreenivasan, Ramshankar and Meneveau1989; de Silva et al. Reference de Silva, Philip, Chauhan, Meneveau and Marusic2013; Krug et al. Reference Krug, Holzner, Marusic and van Reeuwijk2017b), that bears the multiscale properties of turbulence, while its propagation velocity relative to the fluid elements 
$\langle v_{n}\rangle$ is very slow and of the order of the Kolmogorov velocity scale (Holzner & Lüthi Reference Holzner and Lüthi2011). Although the local propagation of the TNTI is of viscous nature, it is well known that the overall entrainment rate is independent of viscosity (Morton, Taylor & Turner Reference Morton, Taylor and Turner1956; Townsend Reference Townsend1966; Tritton Reference Tritton1988; Tsinober Reference Tsinober2009), viz. the Reynolds number. It thus follows that 
$A_{\unicode[STIX]{x1D702}}$ plays a crucial role in setting the entrainment rate, cancelling out the viscosity dependency of 
$\langle v_{n}\rangle$ (Townsend Reference Townsend1966). To date, much of the research on the TNTI and associated entrainment process focused on vorticity transport across the TNTI (Westerweel et al. Reference Westerweel, Fukushima, Pedersen and Hunt2005; Holzner & Lüthi Reference Holzner and Lüthi2011; Silva, Zecchetto & da Silva Reference Silva, Zecchetto and da Silva2018) and little is known about the mechanism that sets the surface area of the TNTI.
In his theoretical work, Phillips (Reference Phillips1972) introduced an equation for the time evolution of the surface area of the TNTI (see § 2), which demonstrated that the growth of the interface area is the result of the sum between a flow stretching term and a curvature/propagation term. Hypothesizing a constant entrainment velocity over the TNTI, he concluded that, on average, the curvature/propagation effect creates TNTI area along the bulges and destroys it in the valleys of the TNTI. While this is an important theoretical finding, the local entrainment velocity is known to vary significantly along the TNTI (Holzner & Lüthi Reference Holzner and Lüthi2011; Wolf et al. Reference Wolf, Lüthi, Holzner, Krug, Kinzelbach and Tsinober2012; Watanabe et al. Reference Watanabe, Sakai, Nagata, Ito and Hayase2014), with a predominance of negative values implying entrainment that alternate with sporadic positive values representing detrainment zones (Wolf et al. Reference Wolf, Lüthi, Holzner, Krug, Kinzelbach and Tsinober2012; Krug et al. Reference Krug, Chung, Philip and Marusic2017a; Mistry, Philip & Dawson Reference Mistry, Philip and Dawson2019). In this work, we evaluate locally the stretching and the curvature/propagation terms on the TNTI in order to assess their role in the time evolution of surface area of the TNTI.
In the last decade, an effort has been made to define the role of coherent structures in the entrainment process. da Silva & dos Reis (Reference da Silva and dos Reis2011) used direct numerical simulation (DNS) data of a turbulent planar jet to show that the large-scale vortices near the TNTI define the shape of the interface area. Related findings were presented by Lee, Sung & Zaki (Reference Lee, Sung and Zaki2017), who used a conditional analysis to show that the surface area of the TNTI increases in the vicinity of large-scale motions of a turbulent boundary layer. More recently, vortical structures near the TNTI have been shown to influence both the intensity of the local entrainment velocity and the mean curvature of the TNTI (Mistry et al. Reference Mistry, Philip and Dawson2019; Neamtu-Halic et al. Reference Neamtu-Halic, Krug, Haller and Holzner2019). This suggests that vortical structures may impact the evolution of the TNTI area. Besides, it was also observed that the coherent vortices near the TNTI distort the mean flow in their proximity (Lee et al. Reference Lee, Sung and Zaki2017; Watanabe et al. Reference Watanabe, Jaulino, Taveira, da Silva, Nagata and Sakai2017), which indicates that they may also influence the stretching of the TNTI. However, to date, the role of coherent flow structures on the time evolution of the TNTI area is largely unknown. Contrary to previous approaches, it is our goal to identify Eulerian vortical structures in a systematic, observer-independent fashion. To this end we detect objective (i.e. frame-independent) Eulerian coherent structures (OECSs) (Haller et al. Reference Haller, Hadjighasem, Farazmand and Huhn2016; Serra & Haller Reference Serra and Haller2016) and elucidate their role on the time evolution of the TNTI area.
Due to their relevance in many geophysical scenarios, turbulent flows with stable stratification have received substantial attention from the scientific community. Examples of such flows include river plumes (MacDonald, Carlson & Goodman Reference MacDonald, Carlson and Goodman2013), cloud-top mixing layers (Mellado Reference Mellado2010) and oceanic overflows (Legg et al. Reference Legg, Briegleb, Chang, Chassignet, Danabasoglu, Ezer, Gordon, Griffies, Hallberg and Jackson2009). For these flows, the entrainment coefficient is known to diminish with increasing ratio between buoyancy and the shear strength of flow, represented by the Richardson number 
$Ri$ (Ellison & Turner Reference Ellison and Turner1959). Recently, it has been demonstrated (Krug et al. Reference Krug, Holzner, Lüthi, Wolf, Kinzelbach and Tsinober2015; van Reeuwijk, Krug & Holzner Reference van Reeuwijk, Krug and Holzner2018) that the reduction of the entrainment coefficient with increasing 
$Ri$ is associated with the decrease of both 
$\langle v_{n}\rangle$ and 
$A_{\unicode[STIX]{x1D702}}$. In particular, Krug et al. (Reference Krug, Holzner, Marusic and van Reeuwijk2017b) showed that the reduction of 
$A_{\unicode[STIX]{x1D702}}$ is caused by the decrease of its fractal scaling exponent, while the scaling range remains largely unaffected. Here, we explore how varying 
$Ri$ affects the role of stretching and curvature/propagation for the time evolution of the TNTI area.
The main scope of the present work is to investigate the mechanisms that continuously produce and destroy the turbulence interface in flows with and without stable stratification with a particular regard to the role played in this process by coherent flow structures and the degree of stratification.
The analysis will be carried out using a DNSs of a temporal wall jet and gravity currents, details of which will be presented in § 2. This is followed by the presentation of the results in § 3, while concluding remarks are given in § 4.
2 Methods
2.1 Direct numerical simulations
In the present paper, we use DNSs of temporal gravity currents and of a temporal wall jet for which there is no stratification. Temporally evolving flows are ideally suited for obtaining converged statistics relatively inexpensively since they are homogeneous in wall-parallel planes (van Reeuwijk et al. Reference van Reeuwijk, Krug and Holzner2018). For the simulations, we solve the Navier–Stokes equations in the Boussinesq approximation with a fourth-order accurate finite volume discretization scheme (Craske & van Reeuwijk Reference Craske and van Reeuwijk2015) on a cuboidal volume of 
$1536\times 1536\times 1152$ cells. Periodic boundary conditions are applied both in the 
$y$ (the spanwise) and the 
$x$ (streamwise) directions. In the 
$z$ (vertical) direction, at the wall (
$z=0$) and at the top of the simulation domain, no-slip respectively free-slip velocity boundary conditions are imposed for the velocity, whereas Neumann (no-flux) boundary conditions are imposed for buoyancy. As schematically represented in figure 1(a), for the initial conditions (indicated by subscript 0), a uniform distribution of both the streamwise velocity 
$u_{0}$ and the buoyancy 
$b_{0}<0$ up to a height 
$h_{0}$ above the bottom wall is used.
Schematic representation of the simulation set-up (a). Time variation of the bulk Richardson number (purple) and the gradient Richardson number (red) (b). Vertical profiles of the mean streamwise velocity (c) and mean buoyancy (d).
The size of the domain in the 
$z$ direction is 
$L_{z}=10h_{0}$, whereas in the 
$x$ and in 
$y$ directions it is 
$L_{x}=L_{y}=20h_{0}$. In order to simulate a sloping bottom, the buoyancy vector 
$\boldsymbol{b}=b\hat{\boldsymbol{g}}$ is tilted at an angle 
$\unicode[STIX]{x1D6FC}$ with respect to the vertical. Here 
$b$ is a scalar with Schmidt number 
$Sc=1$ and 
$\hat{\boldsymbol{g}}=(\sin \unicode[STIX]{x1D6FC},0,-\cos \unicode[STIX]{x1D6FC})g$. In this way, the component 
$b$
$\sin (\unicode[STIX]{x1D6FC})$ drives the flow along 
$x$, while 
$b$
$\cos (\unicode[STIX]{x1D6FC})$ causes a stable stratification in the wall-normal direction. For a more detailed discussion on the DNS concept and numerical configuration we refer to van Reeuwijk et al. (Reference van Reeuwijk, Krug and Holzner2018), whereas the adequacy of the grid resolution is verified in van Reeuwijk, Holzner & Caulfield (Reference van Reeuwijk, Holzner and Caulfield2019).
The different flow cases investigated here differ in the initial Richardson number 
$Ri_{0}=-B_{0}\cos (\unicode[STIX]{x1D6FC})/u_{0}^{2}$, whereas the initial bulk Reynolds number 
$Re_{0}=u_{0}h_{0}/\unicode[STIX]{x1D708}$, where 
$\unicode[STIX]{x1D708}$ is the kinematic viscosity, is kept constant.
To compute the time evolution of 
$Ri=-B_{0}\cos (\unicode[STIX]{x1D6FC})/u_{T}^{2}$, we use the following top-hat definitions:
$$\begin{eqnarray}u_{T}h=\int _{0}^{\infty }\overline{u}\,\text{d}z,\quad u_{T}^{2}h=\int _{0}^{\infty }\overline{u}^{2}\,\text{d}z\quad \text{and}\quad B_{0}=b_{T}h=\int _{0}^{\infty }\overline{b}\,\text{d}z,\end{eqnarray}$$ where 
$B_{0}$ is a conserved quantity in the temporal problem (van Reeuwijk et al. Reference van Reeuwijk, Krug and Holzner2018) and 
$u$ is the streamwise velocity (the overline indicates averaging in wall-parallel planes; the corresponding fluctuations are given by 
$u^{\prime }=u-\overline{u}$). The components of the velocity vector 
$\boldsymbol{u}$ along the 
$y$- and 
$z$-axes are denoted by 
$v$ and 
$w$, respectively. Table 1 summarizes the parameters of the simulations employed in this study. To compute 
$Re_{\unicode[STIX]{x1D706}}$, we average the turbulent kinetic energy 
$e$ and the rate of turbulent dissipation 
$\unicode[STIX]{x1D716}$ in horizontal planes, which were limited at 
$0.3<z/h<1.2$ in order to avoid the influence of the near-wall region (Krug et al. Reference Krug, Holzner, Marusic and van Reeuwijk2017b). Note that the label of the flow cases indicates the value of 
$Ri_{0}$. For the gravity currents (
$Ri11$ and 
$Ri22$), 
$Ri_{0}$ is varied by changing the inclination angle 
$\unicode[STIX]{x1D6FC}$ while keeping the integral forcing 
$\sin (\unicode[STIX]{x1D6FC})B_{0}$ in the 
$x$-direction constant. In addition, we ran a simulation with the buoyancy term switched off, resulting in an unstratified (temporal) wall jet (
$Ri0$) that is driven by initial momentum only. Apart from the § 3.1, where the whole domain was used, results will be based on data over six independent 
$xz$-planes, which are spaced equally in the 
$y$ direction, amounting to 250 snapshots over a period of 
$120h_{0}/u_{0}$. Throughout the paper, the time 
$t$ is normalized by 
$h_{0}/u_{0}$.
Simulation parameters: 
$N_{i}$ and 
$L_{i}$ denote the number of grid points and the size along the 
$i$-direction, respectively. The subscript 0 indicates the inflow parameters. Results for 
$Re_{\unicode[STIX]{x1D706}}$ are averaged over 
$110<t<120$. Here, 
$\unicode[STIX]{x1D707}$ is the kinematic viscosity and 
$\unicode[STIX]{x1D716}$ is the rate of turbulent dissipation.
For a flow characterization, the time evolution of the gradient Richardson number 
$Ri_{g}=(\overline{N}/\overline{S})^{2}$ is shown in figure 1(b). In this definition, 
$\overline{N}^{2}=\text{d}\overline{b}/\text{d}z$ is the buoyancy frequency and 
$\overline{S}=\text{d}\overline{u}/\text{d}z$ is the mean shear. As can be seen from the figure, after an initial transient 
$Ri_{g}$ stabilizes and tends asymptotically towards two different constant values for the gravity current simulations. This behaviour resembles that of the bulk 
$Ri$, however, with slightly lower magnitudes (figure 1b). It is important to note that for 
$Ri_{g}<1/4$, the flow is expected to be shear dominated according to the classification by Mater & Venayagamoorthy (Reference Mater and Venayagamoorthy2014), which has indeed been confirmed for the simulations presented here by Krug et al. (Reference Krug, Holzner, Marusic and van Reeuwijk2017b). Moreover, as can be seen from figure 1(c), when normalized with the top-hat definitions, the mean streamwise velocity profiles of all the flow cases collapse on a single curve. This indicates that, although there are fundamental differences between the stratified and unstratified cases, the structure of the flow is indeed very similar among all the flow cases. It is noteworthy that also the mean buoyancy profile of the gravity currents collapses on a single curve when normalized with the top-hat definitions (figure 1d).
2.2 The TNTI identification and local entrainment velocity
In this study, the position of the TNTI is identified through a threshold on the enstrophy field 
$\unicode[STIX]{x1D714}^{2}=\unicode[STIX]{x1D714}_{i}\unicode[STIX]{x1D714}_{i}$, where 
$\unicode[STIX]{x1D714}_{i}$ is the vorticity vector (Bisset et al. Reference Bisset, Hunt and Rogers2002; Holzner et al. Reference Holzner, Liberzon, Nikitin, Kinzelbach and Tsinober2007, Reference Holzner, Liberzon, Nikitin, Lüthi, Kinzelbach and Tsinober2008; Wolf et al. Reference Wolf, Lüthi, Holzner, Krug, Kinzelbach and Tsinober2012; Silva et al. Reference Silva, Zecchetto and da Silva2018; Neamtu-Halic et al. Reference Neamtu-Halic, Krug, Haller and Holzner2019). Here, the threshold is the same for all flow cases and is set to 
$\unicode[STIX]{x1D714}_{thr}^{2}=10^{-3}{u_{0}}^{2}/{h_{0}}^{2}$, well within the interval of possible values identified by Krug et al. (Reference Krug, Holzner, Marusic and van Reeuwijk2017b) for DNSs with the same code and parameters as the ones presented here. By identifying the TNTI through an isosurface of the enstrophy field, the entrainment velocity can be evaluated using the transport equation for the enstrophy (Holzner & Lüthi Reference Holzner and Lüthi2011; Krug et al. Reference Krug, Holzner, Lüthi, Wolf, Kinzelbach and Tsinober2015). In this case 
$v_{n}$ is given by
$$\begin{eqnarray}v_{n}=-\frac{2\unicode[STIX]{x1D714}_{i}\unicode[STIX]{x1D714}_{j}\unicode[STIX]{x1D61A}_{ij}}{|\unicode[STIX]{x1D735}\unicode[STIX]{x1D714}^{2}|}-\frac{2\unicode[STIX]{x1D708}\unicode[STIX]{x1D714}_{i}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D714}_{i}}{|\unicode[STIX]{x1D735}\unicode[STIX]{x1D714}^{2}|}-\frac{2\unicode[STIX]{x1D716}_{ijk}\unicode[STIX]{x1D714}_{i}\displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{b}_{k}}{\unicode[STIX]{x2202}x_{j}}}{|\unicode[STIX]{x1D735}\unicode[STIX]{x1D714}^{2}|},\end{eqnarray}$$ where 
$\unicode[STIX]{x1D716}_{ijk}$ is the Levi-Civita operator and 
$\unicode[STIX]{x1D61A}_{ij}=1/2(\unicode[STIX]{x2202}u_{i}/\unicode[STIX]{x2202}x_{j}+\unicode[STIX]{x2202}u_{j}/\unicode[STIX]{x2202}x_{i})$ is the strain rate tensor. Throughout the paper, we make use of both fully three-dimensional (3-D) data, as well as two-dimensional (2-D) data from vertical planes. In the case of the 2-D approach, the entrainment velocity is determined through interface tracking, with a procedure similar to the one used by Wolf et al. (Reference Wolf, Lüthi, Holzner, Krug, Kinzelbach and Tsinober2012) in which, 
$v_{n}$ is computed as
$$\begin{eqnarray}v_{n}=v_{I}-u_{n},\end{eqnarray}$$ where 
$v_{I}$ is the local normal velocity of the TNTI and 
$u_{n}=\boldsymbol{u}_{f}\boldsymbol{\cdot }\boldsymbol{n}$, with 
$\boldsymbol{u}_{f}$, the flow velocity at the location of the TNTI and 
$\boldsymbol{n}$, the unit vector normal to the surface of the TNTI pointing towards the turbulent flow region; 
$v_{I}$ is computed by tracking the position of 
$\unicode[STIX]{x1D714}_{thr}^{2}$-isosurfaces in time.
2.3 Equation for the time evolution of the TNTI area
In the present work, we investigate in detail the time evolution of the TNTI area. Initially introduced by Phillips (Reference Phillips1972), the equation for the time evolution of a non-material infinitesimal surface element of area 
$\unicode[STIX]{x1D6FF}A$ reads
$$\begin{eqnarray}\frac{1}{\unicode[STIX]{x1D6FF}A}\frac{\text{d}\unicode[STIX]{x1D6FF}A}{\text{d}t}=(\unicode[STIX]{x1D6FF}_{ij}-n_{i}n_{j})\unicode[STIX]{x1D61A}_{ij}+v_{n}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{n},\end{eqnarray}$$ where 
$\unicode[STIX]{x1D6FF}_{ij}$ is the Kronecker operator, 
$n_{i}$ are the components of 
$\boldsymbol{n}$, the unit vector normal to 
$\unicode[STIX]{x1D6FF}A$, 
$\unicode[STIX]{x1D61A}_{ij}$ is the strain rate tensor and 
$\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{n}$ is the mean curvature of the surface. In this study, 
$\boldsymbol{n}$ points outward towards the non-turbulent side. The first term on the right-hand side of (2.4) is the area stretch term, whereas the second term, hereinafter referred to as the curvature/propagation term, arises from the combined effect of curvature and propagation velocity. Even though (2.4) was initially introduced in the context of studying TNTIs, it has since received more extensive attention from the reactive flow community (Candel & Poinsot Reference Candel and Poinsot1990; Dopazo, Martín & Hierro Reference Dopazo, Martín and Hierro2006). In this field, the equation (2.4) is known as the flame stretch equation.
Since in this work, also 2-D data from vertical planes are employed, we note that by passing from a 3-D to a 2-D approach, the TNTI reduces from a 2-D surface to a one-dimensional (1-D) line and accordingly the symbol 
$A$ is substituted with the symbol 
$l$ for the length of a line element. In this case, the 1-D analogue of (2.4) reads
$$\begin{eqnarray}\frac{1}{\unicode[STIX]{x1D6FF}l}\frac{\text{d}\unicode[STIX]{x1D6FF}l}{\text{d}t}=(\unicode[STIX]{x1D6FF}_{ij}-n_{i}n_{j})\unicode[STIX]{x1D61A}_{ij}+v_{n}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{n},\end{eqnarray}$$ where 
$\boldsymbol{n}$ is the 2-D unit vector normal to the segment 
$\unicode[STIX]{x1D6FF}l$ and 
$\unicode[STIX]{x1D61A}_{ij}$ is the 2-D strain rate tensor.
2.4 Coherent flow structure extraction
For observer-independent vortical structure identification, we employ the recently developed instantaneous vorticity deviation (IVD) technique of Haller et al. (Reference Haller, Hadjighasem, Farazmand and Huhn2016). Derived by Haller (Reference Haller2016) from a new, dynamic version of the classic polar decomposition, the IVD field represents an intrinsic material rotation rate in the fluid. Specifically, the IVD field, defined by the normed deviation of the vorticity vector 
$\unicode[STIX]{x1D74E}(\boldsymbol{x},t)$ from its spatial mean 
$\overline{\unicode[STIX]{x1D74E}}(t)$ over an evolving fluid mass, i.e. by the formula
$$\begin{eqnarray}\text{IVD}(\boldsymbol{x},t)=|\unicode[STIX]{x1D74E}(\boldsymbol{x},t)-\overline{\unicode[STIX]{x1D74E}}(t)|\end{eqnarray}$$provides an observer-independent (objective) local angular velocity for each point of the fluid mass. Outermost tubular surfaces of the IVD, therefore, enable the identification of OECSs in an observer-independent manner, as required for experimentally reproducible coherent structure extraction (Haller Reference Haller2015). To find vortical OECSs in our data set, we use the extraction algorithm developed in Neamtu-Halic et al. (Reference Neamtu-Halic, Krug, Haller and Holzner2019), applied here to the IVD field. In this case, the centre of the vortical structure is represented by a codimension-2 line that is the concatenation of local maxima of the IVD in planes perpendicular to the line itself, whereas the boundary of the structure is the union of the outermost almost-convex isocontours of IVD encircling the local maxima in planes perpendicular to the centreline.
When 2-D data from vertical planes are considered, a 2-D OECS results form the intersection of a 3-D structure with the plane itself. In this case we select only those OECSs with a limited intersection angle with respect to the normal unit vector of the plane. To this end, we impose an upper limit to the ratio between the two eigenvalues of the Hessian of IVD at the location of IVD maxima. The rationale behind this selection is based on the fact that most of the dynamics of tubular vortical structures happens in planes perpendicular to the centreline of the structure.
3 Results
3.1 Time evolution of the TNTI area
Visualization of the stretching term (a–c), curvature/propagation term (d–f) and the time evolution of a non-material infinitesimal element of area 
$\unicode[STIX]{x1D6FF}A$, 
$(1/\unicode[STIX]{x1D6FF}A)\text{d}(\unicode[STIX]{x1D6FF}A)/\text{d}t$ (g–h) of the TNTI for different view angles as captured from a snapshot of the 
$Ri11$ flow case at 
$t=100$. The black arrow indicates the mean flow direction.
A visualization of the TNTI in a sub-domain of the gravity current for the 
$Ri11$ flow case at 
$t=100$ is shown in figure 2. Here, the TNTI is colour coded with the terms of (2.4). In (2.4), positive values of the terms on the right-hand side contribute to the production of surface area of the TNTI, whereas negative values promote its destruction. As can be seen from figure 2(a–f), both 
$(\unicode[STIX]{x1D6FF}_{ij}-n_{i}n_{j})\unicode[STIX]{x1D61A}_{ij}$ and 
$v_{n}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{n}$ act to produce and destroy the TNTI area. In particular, 
$(\unicode[STIX]{x1D6FF}_{ij}-n_{i}n_{j})\unicode[STIX]{x1D61A}_{ij}$ is mainly positive at the leading edges (figure 2b) and negative at the trailing edges of the bulges of the TNTI (figure 2c).
Conversely, 
$v_{n}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{n}$ is particularly active in the valleys of the interface (see e.g. figures 2e and 2f and the zoom of figure 2e), where strong negative values can be observed. This is expected, given that the valleys of the TNTI are regions with high curvature. Moreover, 
$v_{n}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{n}$ appears to have positive values on the bulges, but at a lower intensity as compared to the valleys. The sum of the two terms, which describes the time evolution of a non-material infinitesimal element of area 
$\unicode[STIX]{x1D6FF}A$, is mostly dominated by 
$(\unicode[STIX]{x1D6FF}_{ij}-n_{i}n_{j})\unicode[STIX]{x1D61A}_{ij}$ on the bulges and by 
$v_{n}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{n}$ in the valleys (see figure 2h,i). Interestingly, we note that the sign of the patches, especially for 
$(\unicode[STIX]{x1D6FF}_{ij}-n_{i}n_{j})\unicode[STIX]{x1D61A}_{ij}$ and 
$(1/\unicode[STIX]{x1D6FF}A)\text{d}(\unicode[STIX]{x1D6FF}A)/\text{d}t$, seems to correlate with the geometry of the bulges.
From figure 2, there appears to be a spatial organization of the terms in (2.4) with respect to the TNTI bulges. Neamtu-Halic et al. (Reference Neamtu-Halic, Krug, Haller and Holzner2019) showed that TNTI bulges are populated by OECSs. In figure 3, we show part of the OECSs extracted from a subvolume of the flow field shown in figure 2.
Visualization of the OECSs and the TNTI from 
$Ri11$ at 
$t=100$. The OECSs are represented by blue tubular surfaces (boundaries) surrounding 1-D curves (centres), whereas the TNTI is represented by the yellow open surface. The region above the TNTI is irrotational, whereas below the surface, the flow is turbulent. The OECSs with green boundaries cross the vertical plane at 
$y=9.75$ almost perpendicularly. On the vertical plane at 
$y=9.75$ the IVD field is shown in the red contour plot.
Here, the OECSs are composed of tubular surfaces, that constitute the boundaries of the structures. These surfaces enclose 1-D curves, which represent the centre of the structures. In addition to the OECSs, we also display the nearby TNTI (yellow transparent open surface) along with a vertical 
$xz$-plane (at 
$y=9.75$) colour coded by intensity of the IVD. As observed in Neamtu-Halic et al. (Reference Neamtu-Halic, Krug, Haller and Holzner2019), most of the bulges are filled with OECSs that appear to shape the nearby TNTI. To investigate the local effect of the OECSs on the TNTI area production/destruction, we use the conditional analysis of Neamtu-Halic et al. (Reference Neamtu-Halic, Krug, Haller and Holzner2019), and explore the impact of the coherent structures on 
$(\unicode[STIX]{x1D6FF}_{ij}-n_{i}n_{j})\unicode[STIX]{x1D61A}_{ij}$ and 
$v_{n}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{n}$.
To simplify the analysis and manage the computational cost, we perform the subsequent analysis in two dimensions. Using the selection criterion described in § 2.4, in figure 3, we highlight the 3-D OECSs (green boundaries) that are considered for the further 2-D analysis in the case of the vertical plane at 
$y=9.75$.
Before proceeding, the accuracy of the 2-D approach as compared to the 3-D approach is addressed in terms of probability density functions (PDFs) of the three terms of (2.4). In figure 4, we show the results for the 
$Ri11$ flow case. As can be seen from the figure, in all three PDFs, the two approaches provide very similar results.
In general, the PDF of the stretching term (figure 4a) has a higher positive tail and the overall distribution is slightly shifted towards positive values. Conversely, the PDF of the curvature/propagation term (figure 4b) has a higher negative tail and has a strong peak at 
$v_{n}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{n}=0$. Moreover, the variance of the two aforementioned PDFs is different. In fact, the PDF of the curvature/propagation term presents a much narrower distribution as compared to the PDF of the stretching term. The PDF of the sum of the two terms (figure 4c) shows characteristics of the PDFs of both 
$(\unicode[STIX]{x1D6FF}_{ij}-n_{i}n_{j})\unicode[STIX]{x1D61A}_{ij}$ and 
$v_{n}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{n}$, in that it has a positive peak and a slightly higher negative tail.
The PDFs of stretching (a) and curvature/propagation terms (b) and their sum (c) from 2-D (black) and 3-D data (red) for 
$Ri11$.
3.2 Time evolution of the TNTI area: a 2-D approach
In the following, only results from 2-D data are presented. In figure 5, we show the time evolution of the spatial averages of the terms in (2.4) for each of the flow cases. Note that, 
$(1/\unicode[STIX]{x1D6FF}l)\text{d}(\unicode[STIX]{x1D6FF}l)/\text{d}t$ can be computed locally as a sum of 
$(\unicode[STIX]{x1D6FF}_{ij}-n_{i}n_{j})\unicode[STIX]{x1D61A}_{ij}$ and
$v_{n}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{n}$ only. However, since we consider average values over the entire box, an estimation of the average of 
$(1/\unicode[STIX]{x1D6FF}l)\text{d}(\unicode[STIX]{x1D6FF}l)/\text{d}t$ can be made by taking for 
$\unicode[STIX]{x1D6FF}l$ the entire length of TNTI. This can be formalized according to
$$\begin{eqnarray}\left\langle \frac{1}{\unicode[STIX]{x1D6FF}l}\frac{\text{d}\unicode[STIX]{x1D6FF}l}{\text{d}t}\right\rangle =\frac{1}{\displaystyle \mathop{\sum }_{l}\unicode[STIX]{x1D6FF}l}\frac{\text{d}\left(\displaystyle \mathop{\sum }_{l}\unicode[STIX]{x1D6FF}l\right)}{\text{d}t}=\langle (\unicode[STIX]{x1D6FF}_{ij}-n_{i}n_{j})\unicode[STIX]{x1D61A}_{ij}\rangle +\langle v_{n}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{n}\rangle .\end{eqnarray}$$ As can be seen for all the flow cases, the stretching term is positive on average and decays in time from approximately 
$0.1$ at 
$t=25$ to 
$0.01$ at 
$t=120$. Conversely, the curvature/propagation term is negative on average and its magnitude decays in time similar to the stretching term from approximately 
$-0.1$ at 
$t=25$ to 
$-0.01$ at 
$t=120$. We note that several turbulent time scales were tested for the scaling of these trends, namely, the Kolmogorov time scale 
$(\unicode[STIX]{x1D708}/\unicode[STIX]{x1D716})^{1/2}$, the mean shear time scale 
$\overline{S}^{-1}$, the turbulence time scale 
$e/\unicode[STIX]{x1D716}$ and the integral time scale 
$h/u_{T}$. However, none of these time scales were able to collapse rates in time and across 
$Ri$, thus hinting at a multiscale nature of the terms in (2.4). As 
$Ri$ increases, both the average stretching and curvature/propagation terms decrease. Moreover, the spatial average growth of the TNTI surface, 
$(1/\unicode[STIX]{x1D6FF}l)\text{d}(\unicode[STIX]{x1D6FF}l)/\text{d}t$, results to be very small for all the time steps between 
$t=25$ and 
$t=120$. That is, the two terms on the right-hand side of (2.4) approximately balance each other out for all the time steps shown in figure 5.
Time evolution of the spatial average of stretching term (red), curvature/propagation term (blue) and 
$(1/\unicode[STIX]{x1D6FF}l)\text{d}(\unicode[STIX]{x1D6FF}l)/\text{d}t$ (purple) for 
$Ri0$ (continuous line), 
$Ri11$ (dash-dotted line) and 
$Ri22$ (dashed line).
In order to understand how the two terms on the right-hand side of (2.4) balance out on average, we show in figure 6 the joint PDFs of all possible couples of the terms in (2.4) for 
$Ri0$ (a–c), 
$Ri11$ (d–f) and 
$Ri22$ (g–h).
Joint PDFs of 
$(\unicode[STIX]{x1D6FF}_{ij}-n_{i}n_{j})\unicode[STIX]{x1D61A}_{ij}$, 
$v_{n}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{n}$ and their sum for 
$Ri0$ (a–c), 
$Ri11$ (d–f) and 
$Ri22$ (g–h) for 
$30<t<120$. The corresponding values of the black contour lines increase with logarithmic intervals.
The joint PDF of 
$(\unicode[STIX]{x1D6FF}_{ij}-n_{i}n_{j})\unicode[STIX]{x1D61A}_{ij}$ and 
$v_{n}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{n}$ (figure 6a,d,g) has a vertically elongated shape with a distinguishable horizontally elongated peak. This demonstrates that the stretching term dominates over the curvature/propagation term for small values of 
$(\unicode[STIX]{x1D6FF}_{ij}-n_{i}n_{j})\unicode[STIX]{x1D61A}_{ij}$ and 
$v_{n}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{n}$ (between 
$\pm 0.1$), whereas the curvature/propagation term has longer tails. Furthermore, the joint PDF shows higher probability for positive values of 
$(\unicode[STIX]{x1D6FF}_{ij}-n_{i}n_{j})\unicode[STIX]{x1D61A}_{ij}$ and negative values of 
$v_{n}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{n}$. As 
$Ri$ increases the tails of the curvature/propagation term are reduced, whereas the stretching term shows a slightly broader distribution and a more centred peak. The joint PDF of 
$(\unicode[STIX]{x1D6FF}_{ij}-n_{i}n_{j})\unicode[STIX]{x1D61A}_{ij}$ and 
$(1/\unicode[STIX]{x1D6FF}l)\text{d}(\unicode[STIX]{x1D6FF}l)/\text{d}t$ (sum of the two terms on the right-hand side of (2.4)) is shown in figure 6(b,e,h). A high degree of correlation between the 
$(\unicode[STIX]{x1D6FF}_{ij}-n_{i}n_{j})\unicode[STIX]{x1D61A}_{ij}$ and 
$(1/\unicode[STIX]{x1D6FF}l)\text{d}(\unicode[STIX]{x1D6FF}l)/\text{d}t$ can be noticed for small values (between 
$\pm 0.1$), whereas for higher values the PDF is elongated in the vertical direction, a sign of weaker correlation between the two terms. Again, as 
$R$i increases the tails of 
$(1/\unicode[STIX]{x1D6FF}l)\text{d}(\unicode[STIX]{x1D6FF}l)/\text{d}t$ reduce. The joint PDF between 
$v_{n}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{n}$ and 
$(1/\unicode[STIX]{x1D6FF}l)\text{d}(\unicode[STIX]{x1D6FF}l)/\text{d}t$ is shown in figure 6(c,f,i). In this case, the two quantities correlate very well for intense values, whereas they appear to be almost uncorrelated near the origin. In conclusion, the area growth is mostly driven by the strain term for weak to moderate events that tends to produce interface area. However, the large tails of the curvature/propagation term dominate the extreme events; the term has a negative mean so that, on average, it counterbalances the positive stretching. As 
$Ri$ increases the curvature/propagation term has a narrower distribution, whereas the peak of the stretching term tends to move closer to the origin. A physical interpretation of this latter observation is provided in § 3.3, where we connect the interface evolution to the presence of coherent structures.
3.2.1 Effect of the stable stratification on the production/destruction process of the TNTI area
In the following, we investigate the effect of the stable stratification on the terms of (2.4). Initially, we focus on the stretching term, that in two dimensions can be written as
$$\begin{eqnarray}(\unicode[STIX]{x1D6FF}_{ij}-n_{i}n_{j})\unicode[STIX]{x1D61A}_{ij}=(1-n_{x}^{2})\frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}x}-n_{x}n_{z}\left(\frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}z}+\frac{\unicode[STIX]{x2202}w}{\unicode[STIX]{x2202}x}\right)+(1-n_{z}^{2})\frac{\unicode[STIX]{x2202}w}{\unicode[STIX]{x2202}z}.\end{eqnarray}$$ In figure 7, we show PDFs of the three components of 
$(\unicode[STIX]{x1D6FF}_{ij}-n_{i}n_{j})$ and of the three components of 
$\unicode[STIX]{x1D61A}_{ij}$. As can be seen, a significant effect of the stable stratification can be noticed on the three coefficients 
$(\unicode[STIX]{x1D6FF}_{ij}-n_{i}n_{j})$ of (3.2) (figure 7a–c). While 
$(1-n_{x}^{2})$ shows a higher probability for values close to unity with increasing 
$Ri$ (figure 7a), 
$(1-n_{z}^{2})$ is seen to diminish as the stratification increases (figure 7c). Also, as 
$Ri$ increases, 
$-n_{x}n_{z}$ shows a higher probability for values close to 0. These observations indicate that 
$n_{x}\rightarrow 0$ while 
$n_{z}\rightarrow 1$ as the stratification increases, that is, the interface tends to flatten with increasing 
$Ri$.
The PDFs of the three components of 
$(\unicode[STIX]{x1D6FF}_{ij}-n_{i}n_{j})$ (a–c) for 
$Ri0$ (continuous line), 
$Ri11$ (dash-dotted line) and 
$Ri22$ (dashed line). TNTI surface as obtained from the fractal model of Krug et al. (Reference Krug, Holzner, Marusic and van Reeuwijk2017b) for 
$r=0.3$ (
$Ri0$), 
$r=0.12$ (
$Ri11$) and 
$r=0.09$ (
$Ri22$) (d). The thickness of the lines increases with 
$Ri$. Average value of 
$(1-n_{x}^{2})$ (squares), 
$-n_{x}n_{z}$ (triangles) and 
$(1-n_{x}^{2})$ (circles) from the fractal model (green) and from the DNS data (red) with 
$30<t<120$ (e). The PDFs of the three rate of strain components (f–h) of the stretching term and PDFs of the entrainment velocity (i) and of the curvature of the TNTI (j) for 
$30<t<120$. Average values of 
$(1-n_{x}^{2})\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}x$ (square symbols), 
$-n_{x}n_{z}(\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}z+\unicode[STIX]{x2202}w/\unicode[STIX]{x2202}x)$ (red triangle symbols), 
$(1-n_{z}^{2})\unicode[STIX]{x2202}w/\unicode[STIX]{x2202}z$ (circles) and of 
$\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{n}$ (blue diamonds) and 
$-v_{n}/u_{\unicode[STIX]{x1D702}}$ (blue triangles) against the 
$Ri$ number (j).
To better understand the effect of the stratification on the components of 
$(\unicode[STIX]{x1D6FF}_{ij}-n_{i}n_{j})$ tensor, we use the fractal scaling theory for the geometry of the TNTI (Sreenivasan et al. Reference Sreenivasan, Ramshankar and Meneveau1989). According to the theory, the length of the TNTI depends on 
$l_{i}/l_{o}$, the inner and the outer cutoffs of the scaling range and on 
$\unicode[STIX]{x1D6FD}$, the fractal scaling exponent. In their work, Krug et al. (Reference Krug, Holzner, Marusic and van Reeuwijk2017b) showed that for gravity currents, 
$l_{i}/l_{o}$ is essentially constant for 
$0<Ri<0.22$, whereas 
$\unicode[STIX]{x1D6FD}$ decreases with increasing stratification. Moreover, the authors observed that the convolutions of the TNTI are anisotropic, scaling with 
$l_{sk}$ in the wall-normal direction and with 
$h$ in the streamwise direction. Based on the observation that the ratio 
$r=l_{sk}/h$ decreases with increasing stratification, they implemented a model for 
$\unicode[STIX]{x1D6FD}=f(r)$, that was able to reproduce the trends of the fractal scaling exponent 
$\unicode[STIX]{x1D6FD}$ with increasing 
$Ri$. The model is based on a simple fractal model where in subsequent iterations line segments with length 
$l_{n+1}$ are placed at distance 
$rl_{n}$ from the centre of 
$l_{n}$ (Krug et al. Reference Krug, Holzner, Marusic and van Reeuwijk2017b). Here, we use this model to check quantitatively if the trends observed in figure 7(a–c) can also be related to the anisotropy of the interface bulges. In figure 7(d), we display the geometry of the modelled interface where 
$l_{i}/l_{o}\approx 100$, chosen according to our Reynolds number, and 
$r$ has the values indicated in the caption. As can be noticed, for 
$Ri22$ the fluctuations of the TNTI position in the wall-normal direction are much lower as compared to the wall jet and thus the interface tends to flatten with increasing stratification. A comparison between the fractal model and the DNS data in terms of average value of 
$(\unicode[STIX]{x1D6FF}_{ij}-n_{i}n_{j})$ components against 
$Ri$ is shown in figure 7(e). Notably, the model reproduces rather well the average of 
$(\unicode[STIX]{x1D6FF}_{ij}-n_{i}n_{j})$, especially for 
$(1-n_{x}^{2})$ and 
$(1-n_{z}^{2})$, the two terms that vary more significantly with 
$Ri$. We therefore conclude that, in agreement with the model, the stratification impacts the interface geometry at all ‘active’ length scales of the TNTI between 
$l_{i}$ and 
$l_{o}$.
In figure 7(f–h), we show PDFs of the components of the rate of strain. All three PDFs are only weakly affected by increasing stratification, which is reflected in a moderate increase of the weight of the tail and an associated decrease of the peak at small magnitudes. While in the PDFs of 
$\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}x$ and 
$\unicode[STIX]{x2202}w/\unicode[STIX]{x2202}z$ the skewness reduces (figures 7f and 7h), the PDF of 
$\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}z+\unicode[STIX]{x2202}w/\unicode[STIX]{x2202}x$ is more negatively skewed with increasing 
$Ri$. The latter is consistent with a stronger mean velocity gradient, i.e. smaller 
$h$ at similar 
$u_{T}$, and therefore stronger 
$\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}z$ at increasing 
$Ri$.
The average values of the terms in (3.2) against 
$Ri$ are shown in figure 7(k). Although the PDF of 
$\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}x$ presents a slightly higher negative tail, the average of 
$(1-n_{x}^{2})\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}x$ is positive. This means that there is coupling between high values of 
$(1-n_{x}^{2})$ and positive values of 
$\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}x$. As the stratification increases, 
$(1-n_{x}^{2})\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}x$ increases slightly. The average of 
$-n_{x}n_{z}(\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}z+\unicode[STIX]{x2202}w/\unicode[STIX]{x2202}x)$ is positive as expected from the PDFs in figure 7(b,c). As the stratification increases, 
$-n_{x}n_{z}(\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}z+\unicode[STIX]{x2202}w/\unicode[STIX]{x2202}x)$ increases slightly at 
$Ri=0.11$, to decrease afterwards at 
$Ri=0.22$. In particular, the smaller value 
$-n_{x}n_{z}(\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}z+\unicode[STIX]{x2202}w/\unicode[STIX]{x2202}x)$ at 
$Ri=0.22$ might be related to a higher probability of 
$-n_{x}n_{z}=0$ as compared to the other flow cases, given that the PDF of 
$\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}z+\unicode[STIX]{x2202}w/\unicode[STIX]{x2202}x$ for 
$Ri22$ is comparable to the one of 
$Ri11$ and exhibits even a higher negative tail as compared to the one of 
$Ri0$. Also the average of 
$(1-n_{z}^{2})\unicode[STIX]{x2202}w/\unicode[STIX]{x2202}z$ is positive, meaning that, as for the other terms, there is a coupling of high values of 
$1-n_{z}^{2}$ and positive values of 
$\unicode[STIX]{x2202}w/\unicode[STIX]{x2202}z$. As the 
$Ri$ increases, 
$1-n_{z}^{2}$ tends towards smaller values and the average of 
$(1-n_{z}^{2})\unicode[STIX]{x2202}w/\unicode[STIX]{x2202}z$ decreases.
In conclusion, most of the reduction in the stretching term of (2.4) with increasing stratification is related to a change in the components of 
$(\unicode[STIX]{x1D6FF}_{ij}-n_{i}n_{j})$ tensor as a result of a change in the multiscale geometry of the TNTI which tends to flatten out with increasing 
$Ri$.
Furthermore, in figure 7 we show the impact of the stable stratification on the curvature/propagation term. As can be gleaned from the presented PDFs, the stable stratification reduces the magnitude of both 
$v_{n}$ and 
$\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{n}$. In particular, while the reduction of the 
$v_{n}$ (figure 7i) with increasing 
$Ri$ is well documented (see e.g. Krug et al. Reference Krug, Holzner, Lüthi, Wolf, Kinzelbach and Tsinober2015; van Reeuwijk et al. Reference van Reeuwijk, Holzner and Caulfield2019), it can be seen from figure 7(j) that the stratification reduces also the probability to observe strong convex regions, associated with positive values of 
$\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{n}$. This is not unexpected, given the changes in the geometry of the TNTI discussed above. The average values of the two terms against 
$Ri$ is shown in figure 7(k) and as expected the magnitude of both the terms decrease on average with increasing stratification.
3.2.2 Multiscale nature of the production/destruction of the TNTI area
As observed in figure 2, the positive and the negative patches of the terms in (2.4) appear to correlate with the orientation of the TNTI bulges. Since the TNTI surface has a fractal shape, in figure 8 we investigate the scale dependence of the terms in (2.4). Following the procedure used in Krug et al. (Reference Krug, Holzner, Marusic and van Reeuwijk2017b), we use a box filter to eliminate the effect of the spatial scales smaller than the size of the filter length. To find the position of the interface in the filtered field, we first convert the enstrophy field 
$\unicode[STIX]{x1D714}^{2}$ to a binary field 
$I$, where 
$I=1$ if 
$\unicode[STIX]{x1D714}^{2}>{\unicode[STIX]{x1D714}_{thr}}^{2}$ and 
$I=0$ if 
$\unicode[STIX]{x1D714}^{2}<{\unicode[STIX]{x1D714}_{thr}}^{2}$. We then filter 
$I$ according to 
$\tilde{f}=\int f(\boldsymbol{x}-\boldsymbol{x}^{\prime })G(\boldsymbol{x}^{\prime })\,\text{d}\boldsymbol{x}^{\prime }$, where 
$\tilde{f}$ is the filtered quantity and 
$G$ denotes the kernel of a square box filter of width 
$\unicode[STIX]{x1D6E5}$. We then define the position of the filtered interface as the isocontour 
$I=0.5$. Contrary to filtering the enstrophy field directly, this procedure has the advantage that it preserves the mean position of the TNTI. As highlighted by Krug et al. (Reference Krug, Holzner, Marusic and van Reeuwijk2017b), this is a necessary condition to keep the entrained flux across scales constant. Moreover, to compute the filtered terms in (2.4), we apply the same filter to the streamwise and the wall-normal velocity fields and evaluate the quantities in (2.4) at the location of the filtered interface.
In figure 8, we display the time and space averages of the filtered quantities of (2.4) for different sizes of the filter length 
$\unicode[STIX]{x1D6E5}$, here normalized with the Kolmogorov length scale 
$\unicode[STIX]{x1D702}$. All three flow cases shown here display a similar behaviour, with decaying magnitude of the stretching and the curvature/propagation terms with increasing filter size. Initially for 
$\unicode[STIX]{x1D6E5}/\unicode[STIX]{x1D702}$ smaller than 
${\approx}10$ the decay is very slow. Conversely, for 
$\unicode[STIX]{x1D6E5}/\unicode[STIX]{x1D702}$ larger than 
${\approx}10^{2}$ the magnitude of the stretching and the curvature/propagation terms is negligible, while for 
$\unicode[STIX]{x1D6E5}/\unicode[STIX]{x1D702}$ between 
${\approx}10$ and 
${\approx}10^{2}$ a strong reduction of the two terms can be observed. We note that the limits of the latter region are consistent with the inner (
${\approx}10\unicode[STIX]{x1D702}$) and the outer (
${\approx}0.6{-}0.8h$) cutoffs of the scaling range of 
$A_{\unicode[STIX]{x1D702}}$ observed by Krug et al. (Reference Krug, Holzner, Marusic and van Reeuwijk2017b). Again, as seen in figure 5, when 
$Ri$ increases, the terms on the right-hand side of (2.4) are smaller. The sum of the two terms of (2.4), the variation of the infinitesimal area, remains constant and close to 0 for all the filter sizes. That is to say, the two terms on the right-hand side of (2.4) not only balance out overall, but this balance holds also on a scale by scale basis.
In conclusion, we observe that the difficulty on finding a time scaling for the stretching and the curvature/propagation term discussed in the context of figure 5 is in agreement with the results presented here, in that we showed how these terms are the result of multiscale process. They can therefore not be characterized by a single time scale.
Box and time averages of space-filtered stretching term (red), curvature/propagation term (blue) and the time evolution of a non-material infinitesimal element of area 
$\unicode[STIX]{x1D6FF}l$ (purple) for 
$30<t<120$.
3.3 Impact of OECSs on the production/destruction of TNTI area
A connection between local interface shape and nearby vortical structures was already established in Neamtu-Halic et al. (Reference Neamtu-Halic, Krug, Haller and Holzner2019). This is confirmed in figure 3, where we observed that part of the bulges of the TNTI are filled with OECSs. We now examine how the area change of the TNTI is related to the nearby coherent structures.
As a first step, we focus only on the structures that shape the TNTI, that is, on the structures that are positioned ‘sufficiently’ close to the TNTI. To this end, we compare the size 
$R$ of the structures with the minimum distance 
$\min (d)$ between the boundary of the structure and the TNTI. If the ratio 
$\min (d)/R$ is smaller than a threshold (here fixed at 
$1.5$), the structure is selected and discarded otherwise. The particular value of the threshold was chosen by observing that beyond this limit, the correlation between the structure position and the local interface evolution weakens considerably. In order to determine 
$R$, we fit an ellipse to the boundary of each structure and compute 
$R$ as the mean of the minor and the major axes of the ellipse. To give an impression of the typical size of the structures, we show the PDFs of the size of OECSs near the TNTI, normalized with the Kolmogorov length scale 
$\unicode[STIX]{x1D702}$ for all the flow cases in figure 9(a). As can be seen, all three PDFs have a similar behaviour, increasing rapidly from the minimum value of 
$R/\unicode[STIX]{x1D702}\approx 5$ to the position of the peak at approximately 
$R/\unicode[STIX]{x1D702}\approx 10$. The PDFs then decrease more slightly up to the highest values at about 
$R/\unicode[STIX]{x1D702}\approx 50$. Thus, the size of the OECSs covers a rather broad range of ‘active scales’ over one order of magnitude, which may suggest an important role played by these structures in the multiscale aspects discussed in relation to figure 8.
In order to understand the role of the OECSs on the evolution of the area of the TNTI, we examine whether the relative motion of the coherent structures with respect to the TNTI has an impact on the stretching, and the curvature/propagation term respectively. To this end, for each selected structure, we isolate a region of the TNTI in the neighbourhood of the structure that lies within a distance from the centre of the structure of approximately 
$5R$. This value was chosen in order to select most of the points of the bulge formed by the TNTI in the proximity of the OECS. We then compute the relative velocity 
$\boldsymbol{u}_{f}-\boldsymbol{u}_{c}$ between each of these points of the TNTI and the centre of the OECS. Here 
$\boldsymbol{u}_{f}$ denotes the fluid velocity at the location of the TNTI and 
$\boldsymbol{u}_{c}$ the one at the location of the centre of the OECS. We then project this velocity difference along the connection segment between each point of the TNTI and the centre of the OECS to obtain 
$(\boldsymbol{u}_{f}-\boldsymbol{u}_{c})_{//}$. In figure 9, we show the joint PDF of 
$(\boldsymbol{u}_{f}-\boldsymbol{u}_{c})_{//}$ against the strength of 
$(\unicode[STIX]{x1D6FF}_{ij}-n_{i}n_{j})\unicode[STIX]{x1D61A}_{ij}$ (figure 9b) respectively 
$v_{n}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{n}$ (figure 9c). As can be gleaned from figure 9(b), 
$(\boldsymbol{u}_{f}-\boldsymbol{u}_{c})_{//}$ and 
$(\unicode[STIX]{x1D6FF}_{ij}-n_{i}n_{j})\unicode[STIX]{x1D61A}_{ij}$ are anti-correlated. That is, when 
$(\boldsymbol{u}_{f}-\boldsymbol{u}_{c})_{//}$ is negative, and thus the OECS approaches the TNTI, the stretching term is positive, whereas when 
$(\boldsymbol{u}_{f}-\boldsymbol{u}_{c})_{//}$ is positive the 
$(\unicode[STIX]{x1D6FF}_{ij}-n_{i}n_{j})\unicode[STIX]{x1D61A}_{ij}$ contributes to the TNTI compression. This is not entirely unexpected given that the motion of the OECS towards (away from) the interface implies a compression (stretching) normal to the interface and hence, by incompressibility of water, stretching (expansion) in the tangential plane. A similar behaviour is encountered for the curvature/propagation term and 
$(\boldsymbol{u}_{f}-\boldsymbol{u}_{c})_{//}$ (figure 9c), although less pronounced than in the previous case. Also here, when the OECS approaches the TNTI the curvature/propagation term produces area, whereas when the OECS moves away from the interface, the TNTI area is decreased. Both the joint PDFs shown in figure 9 indicate that the motion of the OECS near the TNTI have an influence on the time evolution of the area of the TNTI.
The PDF of the size (a) of the structures in the proximity of the TNTI. Increasing thickness and transparency of the curves corresponds to increasing 
$Ri$. Joint PDF of the stretching (b) and curvature/propagation terms (c) with respect to the relative velocity between the OECS and the TNTI for 
$Ri11$. The inset in (b) shows schematically where 
$\boldsymbol{u}_{f}$ and 
$\boldsymbol{u}_{c}$ are computed.
3.3.1 Self-similarity of the flow fields around OECSs near the TNTI
In the following, we demonstrate that the flow fields around the OECSs near the TNTI are self-similar with respect to the size of the structures by means of a conditional analysis. In particular, we re-sample the velocity, the enstrophy and the IVD fields onto a uniform grid, frame referenced at the centre of the OECSs, and normalized in 
$x$ and 
$z$ directions with the average size of the structures 
$R$. The rationale here is to have a common frame of reference for all the OECSs and to compare flow fields around OECSs of the same normalized size. By taking the average of the IVD fields around the OECSs, we extract a mean representative OECS, that is, we apply the extraction algorithm to the conditional average of the IVD field. To compute the conditional average of the TNTI position, we apply the same change of the coordinate system to the TNTI position near each selected OECS. By defining a curvilinear coordinate 
$s/R$ along the TNTI which has its origin at the same 
$x$ coordinate of the highest point on the boundary of the OECS, we compute 
$\tilde{x}_{i}(s)$ and 
$\tilde{z}_{i}(s)$, the normalized coordinates of the TNTI position near the 
$i$th OECS. In this way, the mean position of the TNTI is computed by taking the average of 
$\tilde{x}_{i}(s)$ and 
$\tilde{z}_{i}(s)$ over 
$i$ conditioning with respect to 
$s/R$. The same procedure is used also for the terms in (2.4), as well as for the entrainment velocity 
$v_{n}$.
The normalization used here for the spatial coordinates presupposes self-similarity of the velocity and the enstrophy fields in the proximity of the OECSs, as well as of the TNTI position and of the terms in (2.4) with respect to the size of the structures 
$R$. In order to demonstrate this self-similarity, we subdivide the selected structures near the TNTI (of the 
$Ri11$ flow case) in four different groups with similar cardinality based on their size (with the limits indicated in the caption of figure 10) and we conduct a conditional analysis using the spatial normalization indicated above.
In figure 10, we present the results of the conditional analysis. Here, the mean position of the TNTI is colour coded with the conditional average of the time evolution of the surface area of the TNTI. Moreover, in the background the conditional average of the enstrophy field is shown, while the arrows indicate the average velocity field relative to the velocity at the centre of the OECS. Focusing on figure 10(a), the first observation that emerges is that the conditional average position of the TNTI is shaped by the nearby OECSs. Also, the mean flow around the OECS in the frame of reference co-moving with the centre of the structures indicates a rotational fluid motion. The average enstrophy field is particularly high at the centre of the mean OECS, decreasing radially towards the boundaries of the OECS. Furthermore, the conditional average of the infinitesimal area variation (shown in colours at the mean position of the TNTI) displays a clear pattern: positive at the leading edge and negative at the trailing edge. As can be seen in figure 10(b–d), all the groups show identical patterns and there is a striking similarity in all the aspects discussed for figure 10(a). However, as the size of the structures increases (from figure 10a to figure 10d), both the intensities of the enstrophy field inside the boundaries of the OECS and of the flow field around the OECS decrease. This means that the smaller the structures are, the faster is their rotational motion. Moreover, also the typical pattern shown by the conditional average of the time evolution of the TNTI appears to weaken, when the size of the structures is larger. In conclusion, the results shown in figure 10 demonstrate clearly that the vortical structures near the TNTI exhibit a scale-invariant behaviour.
Mean OECS and conditional average position of the TNTI for 
$Ri11$. The boundary of the mean OECS is constituted by the red closed curve. The open curve is the conditional average position of the TNTI, colour coded with the time evolution of a non-material infinitesimal element of area 
$\unicode[STIX]{x1D6FF}l$. The direction and the size of the vectors represent the conditional average velocity field relative to the velocity at the centre of the OECSs, whereas in the background the conditional average of the enstrophy field is shown. The conditional analysis is limited to structures near the TNTI with respectively 
$R/\unicode[STIX]{x1D702}<10.5$ (a), 
$10.5<R/\unicode[STIX]{x1D702}<12.75$ (b), 
$12.75<R/\unicode[STIX]{x1D702}<15$ (c) and 
$R/\unicode[STIX]{x1D702}>15$ (d).
Same conditional average shown in figure 10 for 
$Ri11$. In this case, the conditional average position of the TNTI is colour coded with the average entrainment velocity (a), stretching (b) and curvature/propagation terms (c) and (d) the time evolution of a non-material infinitesimal element of area 
$\unicode[STIX]{x1D6FF}l$.
3.3.2 Conditional analysis of the time evolution of the TNTI
The effect of the OECSs on the terms of (2.4) is investigated in the following through the conditional analysis presented in § 3.3.1. In figure 11, we present the results of the conditional analysis for 
$Ri11$ in which all the selected structures near the TNTI are considered. In particular, in figure 11(a), the mean TNTI is colour coded with the conditional average of the entrainment velocity conditioned with respect to the position on the TNTI. As can be seen, 
$v_{n}$ varies along the TNTI and in particular, it is negative at the leading edge, where it reaches values as high as 
$v_{n}/u_{\unicode[STIX]{x1D702}}=-1$, whereas it decreases (in magnitude) at the trailing edge. This result is consistent with the findings from the experimental work by Neamtu-Halic et al. (Reference Neamtu-Halic, Krug, Haller and Holzner2019), where a very similar pattern for 
$v_{n}$ was observed. In figure 11(b), we display the conditional average of the stretching term of (2.4). Here, the mean stretching term has a clear pattern along the TNTI. Indeed, it has positive values at the leading edge respectively negative values at the trailing edge. On the tails, far from the centre of the structure, it shows smaller and positive values recovering the unconditioned average value. The surface of the TNTI shown in figure 11(c) is color-coded the conditional average of the curvature/propagation term of (2.4). The mean curvature/propagation term exhibits a similar behaviour as the stretching term, being slightly positive at the leading edge and negative at the trailing edge. The sum of the two terms is plotted in figure 11(d). The average infinitesimal area growth rate shown here is positive at the leading edge and negative at the trailing edge. That is, the interface area is produced at the leading edge and destroyed at the trailing edge, whereas far from the structure is negligible reaching the unconditioned average value.
Conditional average of the stretching (blue) and the curvature/propagation term (red) and of their sum (purple) in the proximity of the OECS against the curvilinear abscissa 
$s/R$ for 
$Ri0$ (continuous line), 
$Ri11$ (dash-dotted line) and 
$Ri22$ (dashed line). In the inset, a schematic of the conditional average flow direction around the OECS in the proximity of the interface is shown. The arrows indicate the flow direction between the boundary of the OECS and the TNTI (blue) and in the turbulent region below the OECS (black). ‘
$C$’ denotes the compression, whereas ‘
$S$’ denotes the stretching.
In the following, a comparison of the conditional averages of terms in (2.4) among the different flow cases is presented. For an easier quantitative analysis, in figure 12 the conditional averages are shown against the curvilinear abscissa 
$s/R$. As can be noticed, all three terms of (2.4) have similar trends independently of the flow cases. In particular, moving from 
$s/R=-6$ in the direction of increasing 
$s/R$, the stretching term decreases, reaching a minimum at about 
$s/R\approx -2$. Continuing in the direction of increasing 
$s/R$, 
$(\unicode[STIX]{x1D6FF}_{ij}-n_{i}n_{j})\unicode[STIX]{x1D61A}_{ij}$ starts to increase and after an inflection point at 
$s/R=0$, it reaches a maximum at 
$s/R\approx +2$ and then it decreases again. The behaviour just described is more intense for the stratified cases, where the minimum of 
$(\unicode[STIX]{x1D6FF}_{ij}-n_{i}n_{j})\unicode[STIX]{x1D61A}_{ij}$ at 
$s/R\approx -2$ is negative. For 
$Ri0$, the 
$(\unicode[STIX]{x1D6FF}_{ij}-n_{i}n_{j})\unicode[STIX]{x1D61A}_{ij}$ has the same behaviour but it never reaches negative values. Interestingly, the inflection point of 
$(\unicode[STIX]{x1D6FF}_{ij}-n_{i}n_{j})\unicode[STIX]{x1D61A}_{ij}$ happens exactly at 
$s/R=0$. A possible mechanism that generates the tangential stretching/compression at the leading/trailing edge of the TNTI surface near the vortical structures is schematically depicted in the inset of figure 12. Here the arrows indicate the conditional average flow direction of the fluid motion between boundary of the OECS and the TNTI (blue) and in the turbulent region below the OECS (black). As can be seen, at the leading edge the flow tends to stretch fluid elements, when projected on the tangent to the TNTI surface. Conversely, at the trailing edge, fluid elements are compressed along the direction tangent to the TNTI. In the mechanism just described, the faster the fluid layers beneath the structures, the higher the difference between the relative maximum and minimum of the stretching term. As 
$Ri$ increases, faster and more horizontal layers of fluid are known to form in the turbulent region. These layers might be responsible for the observed difference among the flow cases. We also note that in the conditional average stretching shown in figure 12, the positive contribution of the larger scales is superimposed to the mechanism just described. This is why, for 
$Ri0$ the relative minimum of the stretching term never reaches negative values. Finally, it is worth noting that the formation of the positive–negative double peak seen for 
$Ri22$ is compatible with the more symmetric stretching observed in figure 6.
Starting from 
$s/R=-6$ where it is negative and moving into direction of positive 
$s/R$, the conditional average of the curvature/propagation term increases, reaching a positive maximum at 
$s/R\approx 1$, and decreases, turning again negative at about 
$s/R=2$. To explain this behaviour, we notice that on both sides of a bulge the curvature is positive (valleys), whereas the entrainment velocity is negative (ambient fluid is entrained). Thus their product is negative and it contributes to the destruction of the area of the TNTI. This is why, negative values of the curvature/propagation term are observed at both sides of the OECSs. Moving from downstream to the leading edge, the curvature of the interface changes the sign (at about 
$s/R\approx 2$) becoming negative (concave shape as seen from the turbulent side). The entrainment velocity is still negative here, and thus the curvature/propagation term becomes positive. A weak effect of the stratification can be noticed also here where a higher positive value for the maximum of 
$v_{n}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{n}$ can be observed for the unstratified case.
The sum of the two terms, the time variation of the infinitesimal area, shows a profile that is similar to the stretching term. Indeed in the proximity of the structure, 
$(\unicode[STIX]{x1D6FF}_{ij}-n_{i}n_{j})\unicode[STIX]{x1D61A}_{ij}$ dominates the time evolution of the interface area, dominating over 
$v_{n}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{n}$. Thus, the time variation of the infinitesimal area is positive at the leading edge, it has an inflection point at about 
$s/R\approx 0$ and it becomes negative at the trailing edge. That is, the area of the TNTI is produced at the leading edges and destroyed at the trailing edges. This behaviour is clearly consistent with the trends observed qualitatively in figure 2.
4 Concluding remarks
In this paper, we investigated the production/destruction process of the TNTI area of gravity current flows. We showed that a curvature/propagation term, which is particularly active in the valleys of the TNTI (figure 2), contributes on average to the destruction of the TNTI area (figure 5), while this is counterbalanced by a flow stretching which on average produces interface area (figure 5). In particular, the stretching drives the time evolution of the TNTI area for weak to moderate events, whereas the curvature/propagation term dominates the extreme events (figure 6). Very similar results have been found by Wang et al. (Reference Wang, Hawkes, Chen, Zhou, Li and Aldén2017) in the case of a premixed jet flame. In their work, the authors report that the tangential strain produces interface in regions with low curvature, whereas the curvature/propagation term destroys flame surface area in high curvature regions.
The multiscale aspects of the production/destruction of the TNTI area were investigated by filtering the scales smaller than a filter length 
$\unicode[STIX]{x1D6E5}$ (figure 8). We showed that the stretching and the curvature/propagation terms balance each other at all scales and that the magnitude of the two terms decreases with increasing 
$\unicode[STIX]{x1D6E5}$. Recently, Mistry et al. (Reference Mistry, Philip, Dawson and Marusic2016) showed that the mass-flux rate of entrained fluid across the TNTI is constant across all the scales 
$\unicode[STIX]{x1D6E5}$, as envisaged by Meneveau & Sreenivasan (Reference Meneveau and Sreenivasan1990). That is, the reduction of 
$A_{\unicode[STIX]{x1D702}}$ with increasing 
$\unicode[STIX]{x1D6E5}$ is balanced by the enhancement of 
$v_{n}$ at the larger scales. In our case, 
$v_{n}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{n}$ reduces with increasing filter length, which means that the increment of 
$v_{n}$ is smaller than the rapid decrease of curvature with increasing 
$\unicode[STIX]{x1D6E5}$.
We showed that increasing stratification reduces the averages of both the stretching and the curvature/propagation terms (figure 5), while maintaining the same trends for their time variation (figure 5) and their filtered values (figure 8). In particular, we showed that the reduction of the stretching term is largely attributable to changes in 
$\unicode[STIX]{x1D6FF}_{ij}-n_{i}n_{j}$ tensor (figure 7) and we associated this with the tendency of the interface to flatten with increasing stratification. Indeed, as shown by Krug et al. (Reference Krug, Holzner, Marusic and van Reeuwijk2017b), for the same flows investigated here, the convolutions of the TNTI scale with the shear length scale 
$l_{sk}=e^{1/2}/\overline{S}$ in the vertical direction respectively with 
$h$ in the streamwise direction. While the range of the length scales of the TNTI convolutions impacted by the stable stratification remains the same, it is the growing anisotropy implied by different scalings that modify the geometry of the interface with increasing 
$Ri$. Moreover, this change of the TNTI geometry is also compatible the observation that the positive tail of the interface curvature reduces with increasing 
$Ri$, which together with the entrainment velocity reduction explains the reduction of the curvature/propagation term.
The local effect of the coherent structures on the TNTI area production/destruction process was investigated through a conditional analysis similar to Neamtu-Halic et al. (Reference Neamtu-Halic, Krug, Haller and Holzner2019). The conditional analysis showed that both stretching and the curvature/propagation terms produce TNTI area at the leading edge and to destroy it at the trailing edge of the TNTI in proximity of vortical structures. In particular, we inferred that the behaviour of the stretching term is related to the mean rotational motion induced by the vortical structures (figure 12). A similar mechanism for the tangential stretch of a premixed flame was described in Sinibaldi et al. (Reference Sinibaldi, Driscoll, Mueller, Donbar and Carter2003), where the authors observed that toroidal vortices near the flame boundaries induce a rotational motion that stretches the boundaries on one side of the vortices and compresses it on the other.
Using the conditional analysis, we demonstrated the existence of a scale invariant behaviour of the vortices near the TNTI (figure 10). A self-similarity of the shapes of the TNTI has been observed a long time ago (Sreenivasan et al. Reference Sreenivasan, Ramshankar and Meneveau1989); however, a similar observation for the coherent structures near the TNTI is missing up to date in the literature. Moreover, we showed that also the conditional average of the time evolution of the TNTI is self-similar with respect to the size of the OECS. In particular, the time evolution of the TNTI area exhibits lower values for increasing size of the OECS. In figure 12, we displayed that near the vortical structures, the time evolution of the interface is mostly dominated by the stretching term. This together with the observation that the rotational motion of the OECS is slower for increasing size of the OECSs (figure 10) lends support to the interpretation that the stretching mechanism of the interface is related to the rotation of the OECSs.
In conclusion, to our knowledge, the detailed analysis presented here constitutes a first tentative to describe the evolution of the surface area of the TNTI in the case of turbulent stably stratified shear flows. This work may motivate further studies into the production/destruction mechanisms of the TNTI of other types of turbulent shear flow, such as turbulent jets and wakes and turbulent boundary layers. In particular, it may be of interest to understand whether the vortical structures near the TNTI of these flows have a similar impact on the time evolution of the TNTI and universally apply to all shear flows with an interface. A further issue that may deserve further study is how stratification acts to suppress stretching at large inertial scales.
Acknowledgements
We are grateful for financial support from DFG priority programme SPP 1881 under grant number HA 7497/1-1. J.-P.M. and M.v.R. were supported by the EPSRC project Multi-scale Dynamics at the Turbulent/Non-turbulent Interface of Jets and Plumes (grant number EP/R043175/1) and the UK Turbulence consortium (grant number EP/R029326/1).
Declaration of interests
The authors report no conflict of interest.
