Criteria of tracking vortex surfaces in turbulent-like flows

We propose criteria of tracking vortex surfaces in complex flows based on the vortex-surface field (VSF). The criteria characterize the accuracy and Lagrangian tracking performance of the numerical VSF solution, and determine the time period when the vortex surface tracking is satisfactory. Moreover, we develop a turbulent-like flow combining large-scale coherent structures in the Taylor–Green flow and small-scale turbulent structures in homogeneous isotropic turbulence (HIT). From tracking of vortex surfaces during the effective tracking period, we find that the imposed HIT disturbance significantly wrinkles vortex surfaces. Subsequently, the wrinkled vortex tube with large vorticity magnitude tends to be further twisted, contributing to energy cascade, while the wrinkling is mitigated in the region with small vorticity magnitude.


Introduction
The vortex-surface field (VSF) is a Lagrangian-based method to identify and characterize coherent vortical structures [1,2]. The VSF isosurface is a vortex surface consisting of vortex lines. From the extended Helmholtz theorem [3], the VSF isosurfaces of the same threshold at different times can have strong coherence, facilitating the tracking of vortex surfaces.
As a general flow diagnostic tool, the numerical VSF solution is calculated from a given flow field by solving a pseudo-transport equation driven by the frozen, instantaneous vorticity [4,5]. From an initial VSF, the two-time method [6] is applied to calculate the temporal evolution of VSFs. By post-processing of a time series of numerical solutions, the VSF can elucidate mechanisms in various aerodynamic flows with essential vortex dynamics, such as turbulence [7] and transition [8]. For example, the VSF reveals the continuous rolling-up and twisting of vortex tubes in Taylor-Green (TG) flows [6], and the complex network of tangling vortex tubes in homogeneous isotropic turbulence (HIT) [5], instead of the visual "breakdown" of worm-like structures identified by Eulerian vortex criteria [9,10].
Although the evolution of vortex surfaces has been reported in relatively simple transitional flows, such as the TG flow [6,11] and canonical wall-bounded transitional flows [12,13], the tracking of vortex surfaces in fully developed turbulence, such as HIT, remains an open problem. As illustrated in Fig. 15 in Ref. [5], the successful vortex surface tracking can be valuable to study the physical picture of energy cascade [7,14,15] and vortex-based models [16,17] in turbulence.
The challenges of tracking vortex surfaces in turbulent-like complex flows are rooted in the breakdown of the Helmholtz theorem in real flows, which can cause the non-existence and non-uniqueness of VSF solutions. The former issue has been partially resolved using the two-time method. The error of numerical VSF solutions can be controlled under 5% in various applications. On the other hand, the latter issue can hinder an effective tracking of vortex surfaces [3] due to the correction substep in the two-time method, and there lacks a criterion to quantify this issue. In particular, the numerical dissipation in VSF calculation has two sides. It can regularize the numerical solution to avoid the formation of nearly-singular structures, while it can also artificially distort vortex surfaces to reduce time coherence in the surface tracking.
In the present study, we address the issue on the comprehensive assessment of vortex surface tracking by proposing new criteria in the VSF calculation, and develop a turbulent-like flow combing the large-scale coherent structures and small-scale turbulent structures to assess VSF solutions with moderate computational resources. The tracking criteria determine the effective tracking period. The evolution of VSF isosurfaces during this period demonstrates timewise evolution process of coherence structures, and reveals physically interesting vortex dynamics on the mechanism of transition and energy cascade.
The outline of this paper is as follows. In Section 2, we present the numerical method for the direct numerical simulation (DNS) and the construction of the turbulent-like flow. In Section 3, we review the theory and numerical method of VSF calculation. In Section 4, we propose criteria of tracking vortex surfaces based on the numerical VSF solution. The criteria are applied to the VSF evolution in the turbulent-like flow to investigate vortex dynamics in scale cascade in Section 5. Some conclusions are drawn in Section 6.

DNS
The vortex surface tracking is applicable to general flows [6,11]. Without loss of generality, we restrict our discussion to constant-density, incompressible viscous flows in the present study. The velocity field u(x, t) is governed by the Navier-Stokes (NS) equations where x denotes spatial coordinates, t the time, ρ the density, p the pressure, and ν the kinematic viscosity. We carried out the DNS with the standard pseudo-spectral method [18,19]. The spatial domain of the flow field is a cube of side L = 2π, and it is discretized on uniform grid points N 3 . The two-thirds truncation method with the maximum wavenumber k max ≈ N/3 is applied to remove aliasing errors. The Adams-Bashforth method is used for time advancement. The time step satisfies that the Courant-Friedrichs-Lewy (CFL) number is less than 0.2 to ensure the numerical stability and accuracy.

Construction of the turbulent-like flow
We propose a turbulent-like flow in which the vortex tracking criteria can be effectively assessed. This turbulent-like flow is a TG flow imposed by a HIT perturbation before transition. It consists of deterministic large-scale structures of the TG flow and stochastic small-scale structures of HIT.
The initial velocity is where u TG and u HIT denote the TG and HIT velocity fields, respectively, and η M ≥ 0 denotes the strength of the HIT perturbation. Both u TG and u HIT are calculated by independent DNS cases with the time t TG for the TG flow and t HIT for the HIT starting from their own initial velocity fields u TG (t TG = 0) and u HIT (t HIT = 0), respectively. The hybrid velocity in Eq. (3) consists of u TG at a particular time t * TG and u HIT at another particular time t * HIT . For the DNS of the TG flow, we set u TG (t TG = 0) = (sinxcosycosz, −cosxsinycosz, 0) [20,21], and add the HIT disturbance at t * TG = 4.5 around the beginning of TG transition when the vortex reconnection occurs [6]. The Reynolds number is set as Re ≡ 1/ν = 400. The imposed HIT perturbation can trigger the emergence of convoluted vortex surfaces in the transition.
For the DNS of HIT, u HIT evolves from an initial random velocity field u HIT (t HIT = 0) satisfying a prescribed energy spectrum E(k, t = 0) ∼ k 4 e −(k/k p ) 2 [22,23] with the wavenumber magnitude k and the peak wavenumber k p = 4 [5]. The initial total energy E 0 ≡ |u HIT (t HIT = 0)| 2 /2 is set to be unity, where · denotes the volume average over the computational domain . We choose t * HIT = 1.3 when the HIT becomes fully developed with the maximum mean dissipation rate in the temporal evolution and the Taylor-scale Reynolds number 54 [5]. The disturbance strength is chosen as η M = 0.05 and 0.1 so that the flow is featured by both large-scale evolving TG structures and small-scale turbulent-like structures.

DNS statistics
We set up three cases with the same Re = 400 and three η M = 0, 0.05, 0.1 listed in Table 1. Note that case M0 with η M = 0 is equivalent to the regular TG flow. Figure 1 plots the temporal evolution of the total energy E tot ≡ E(k)dk, and the mean dissipation rate ε ≡ 2νk 2 E(k)dk. We observe that E 0 increases with η M and E tot decays with time in all cases. Meanwhile, ε first generally grows, then reaches a plateau around t = 1.5, and gradually decays at late times. The growth of ε implies a transition process with the vortex reconnection and the emergence of small-scale vortical structures [5,6]. The notable dip for case M2 in the early stage appears to be due to a transient process for reconciling the TG base flow and the imposed HIT disturbance.  Figure 2 plots the energy spectra for u(t = 0) and u HIT (t HIT = t * HIT ). The imposed HIT disturbance enhances E(k) in cases M1 and M2. In particular, the small-scale random turbulence structures smooth E(k) at moderate and large wavenumbers. The temporal evolution of E(k) in cases M1 and M2 in Fig. 3 indicates the emergence of turbulent-like structures at moderate wavenumbers, and E(k) in M2 is slightly larger than that in M1 at moderate k. At t = 2 when ε peaks, E(k) at a small range of k satisfies the −5/3 scaling for low Re = 400, and this range is expected to be broadened with increasing Re.

VSF theory
The VSF φ v , a global smooth scalar field, is defined to satisfy the constraint where ω ≡ ∇ × u denotes the vorticity. The uniqueness and existence of the solution to Eq. (4) have been extensively discussed [1]. An exact VSF solution can be obtained for a flow field with vanishing helicity density [24], and an approximate VSF solution can be calculated from an arbitrary flow field [5]. From the Lagrangian view, the surface tracking in the narrow sense is described by the flow map of a specified surface. By the Helmholtz and Ertel theorems, the tracking of vortex surfaces is well defined in an ideal flow which is inviscid, barotropic, and with conservative body forces. It is equivalent to track a specified isosurface of the VSF governed by the Lagrangian advection equation On the other hand, the vortex surface tracking is generally ill-defined in non-ideal flows due to the breakdown of the Helmholtz theorem. One attempt to resolve this issue introduced a virtual velocity [3], instead of the physical one, to advect the VSF in Eq. (5), but it is only applicable to very special flows. Therefore, a precise definition of tracking of a vortex surface is still an open theoretical problem, but this issue can be reasonably approximated and regularized using numerical methods for practical interests.

VSF calculation
In non-ideal flows, the temporal evolution of the VSF can be numerically solved by the two-time method [6]. This numerical regularization method has been demonstrated to be effective to evolve the VSF and reveal the Lagrangian-like evolution of vortex surfaces in a range of non-ideal flows, e.g., incompressible and compressible viscous TG flows [6,11], reacting flows [25], and magnetohydrodynamic (MHD) flows [3].
Each time step of the two-time method consists of two substeps. In the prediction substep, the temporal VSF φ * v is advected by the physical velocity as a Lagrangian scalar as Then in the correction substep, φ * v is transported by the frozen vorticity with a pseudo time τ as where T τ is the maximum value of pseudo-time. The theoretical analysis and numerical experiment showed that φ * v converges to an approximate VSF solution at large T τ [6]. The initial condition φ v (x, t; τ = 0) for Eq. (7) is obtained from the result in the prediction substep. At the end of each time step, (6) and (7) have no diffusion terms, so in principle nearly singular structures can form in φ * v at large times or pseudo times. A dissipative numerical scheme such as the weighted essentially non-oscillatory (WENO) scheme [26,27] was used for regularization. However, the smoothing of the nearly singular structures can also make φ * v deviate from an accurate one. Thus, we need to carefully choose T τ to balance the convergence and the numerical dissipation of φ * v in Eq. (7), and we typically set 10 t ≤ T τ ≤ 50 t. In the numerical implementation of solving Eqs. (6) and (7), we use the second-order TVD Runge-Kutta method for time advancement and the fifth-order WENO for spatial discretization of the convection terms. Moreover, the choice of spatial discretization schemes is discussed in Appendix A. The time and pseudo time steps are selected by the CFL numbers based on the velocity and vorticity magnitudes [6] less than 0.2, respectively. In general, the computational cost of VSF calculation is higher than that for solving the NS equations, and it depends on the desired solution quality, maximum value of vorticity magnitude, and complexity of vortex lines in flow fields [4].
In case M0, the initial VSF φ 0 at t = 0 is obtained as φ 0 = φ TG (t TG = t * TG ) from the evolution of the VSF φ TG in the TG flow with the exact VSF φ TG (t TG = 0) = (cos 2x − cos 2y) cos z. In cases M1 and M2, φ 0 is calculated using the pseudo-transport in Eq. (7) with the local optimization [5] and T τ = 2400 t. The initial guess of Eq. (7) is φ TG (t TG = t * TG ), and the local optimization is applied to each pseudo-time step. This hybrid method balances the accuracy and smoothness of the converged VSF solution.

VSF deviation
As discussed in Section 3.1, the tracking of vortex surfaces is perfect using Eq. (5) in ideal flows, whereas the assessment of the tracking becomes complex in non-ideal flows. On the accuracy of the VSF solution, Yang and Pullin [6] proposed a metric to quantify the averaged deviation of a numerical VSF solution from an exact one with the cosine of an angle between the VSF gradient and the vorticity, characterizes the local VSF deviation. Former VSF studies [5,6] suggested that a satisfactory VSF solution should have φ ≤ 5% ∼ 10%.
On the other hand, it is not sufficient to conclude that a temporally evolving VSF solution with a small φ must be a good candidate for tracking vortex surfaces due to the nonuniqueness of the VSF solution [1]. An extreme example provided in a MHD flow (see in Fig. 8 in Ref. [3]) shows that the vortex surface should be time invariant, but the VSF isosurface evolves among nonunique VSF solutions with φ = 0.
We sketch this issue in Fig. 4. Two typical vortex surfaces represented by circular dashed lines at a time t evolve into the two represented by elliptical dashed lines at a later time t + t. An ideal vortex surface tracking should show the one-to-one map between the two inner loops or the two outer loops at the two times in a Lagrangian manner. There are two types of discrepancies in the tracking. For the inner vortex surface at t (thick solid line), two isosurfaces of numerical VSF solutions are sketched at t + t. The upper one is attached across the inner and outer vortex surfaces with satisfactory φ and poor Lagrangian tracking performance, i.e., the solid loop obviously deviates from the inner dashed loop. On the contrary, the bottom one is oscillating around the inner vortex surface with large φ and satisfactory Lagrangian tracking performance. Thus, besides φ , we need another metric to quantify the Lagrangian tracking performance in the overall assessment of the VSF solution.

Volume overlap ratio of VSF isosurfaces
The evolution of vortex surfaces is characterized by the VSF isosurfaces of a particular contour value φ v = ϕ. Although ϕ is expected to be constant, it is adjusted at different Fig. 4 Sketch of different surface tracking performances from a given time t (left) to a later time t + t (right). The blue dashed loops denote two vortex surfaces in the ideal tracking, and red solid loops denote a VSF isosurface from a numerical VSF solution. The upper right solution has satisfactory accuracy but poor Lagrangian tracking performance, and the bottom right one has poor accuracy but satisfactory tracking performance times due to the numerical dissipation. In this implementation, ϕ is implicitly determined from [11] by searching the isocontour values φ v = ϕ at a given time t and φ v = C φ at the initial time which correspond to the same fluid volume Here, the fluid volume within an isosurface of φ v = ϕ is calculated by the integration over the entire computational domain as where is the discretized Heaviside function [28] with a smoothing parameter δ = 1 × 10 −6 . The previous VSF visualization shows that the isosurface of φ v = ϕ can reasonably describe the continuous temporal evolution of a particular vortex surface, and this VSF isocontour selection is consistent with the volume conservation implied by the Helmholtz theorem in the limit of in ideal flows [11].
To quantify the tracking performance from numerical VSF solutions, we propose a volume overlap ratio where the volume overlap is calculated by Here, φ v is the VSF solution being examined, and φ ref v denotes an exact reference VSF solution. In the implementation, φ ref v is obtained as a numerical VSF solution with high grid resolution, so it is presumed to be accurate enough with good convergence and low numerical dissipation. The isocontour values of ϕ and ϕ ref satisfy Thus, R V characterizes the normalized volume overlap ratio between the volumes within the isosurfaces of the candidate VSF solution φ v = ϕ and the reference VSF solution For example, Fig. 5 plots cross-sections of VSF solutions on the y-z slices at x = 3 and t TG = 5.0 in the TG flow for Re = 400 with N = 128 (red for the VSF solution to be tested) and N = 512 (blue for the reference VSF solution) along with their overlap region (shaded), where Fig. 5c combines Fig. 5a and b. The isocontour value corresponds to V /8 and the volume overlap ratio is R V = 0.44, where V denotes the volume of the whole computational domain.
The volume overlap ratio characterizes the time coherence of the vortex surface tracking based on VSF solutions. The best tracking performance has R V = 1. To balance reasonably good time coherence and moderate computational cost, we propose a criterion R V ≥ 0.5 for satisfactory performance of vortex surface tracking. Note that the VSF solution is closer to the reference one while the VSF grid resolution is more demanding for larger R V .

Conservation of numerical VSF solutions
In practice, R V may be not available due to the lack of the reference VSF solution with high grid resolution, so we need to find a surrogate metric to assess the Lagrangian tracking performance of numerical VSF solutions. In the present VSF evolution, the VSF isosurface is volume-preserving under the selection of isocontour levels in Eq. (10), but it can be artificially distorted by the dissipative numerical regularization.
The distortion is related to the numerical dissipation, which can be quantified by Note that R φ = 1 for the ideal vortex surface tracking based on Eq. (5). In similar front tracking problems, such as in the level set method [29,30], the conservation degree of the numerical solution of the tracer scalar governed by the equation in the form of Eq. (5) is an important metric to assess the tracking performance. Thus, we examine the correlation of R φ and R V to demonstrate that whether R φ is suitable to characterize the vortex surface tracking.
In the VSF calculations for the TG flow with Re = 400, we use the VSF solution with N = 512 as the reference one, and collect R V and R φ at different times and with different grid resolutions. Note that in order to avoid the grid resolution effect of velocity field, we set the same number of grid points 512 3 for DNS but different N for VSF calculations. In Fig. 6, the data points on the same curve from right to left are obtained at t TG = 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0. In this time period, the TG flow has a transition to the turbulent-like state. The isocontour values correspond to V φ = V /16 and V /4.
We find that R φ and R V are positively correlated, indicating the effectiveness of R φ for characterizing the surface tracking performance. Moreover, the correlation depends on the grid resolution and the isocontour level. It is noted that the volume V φ within the isosurface of φ v = C φ decreases with the increasing C φ . For the same R φ , R V increases with N and V φ . This suggests that the isosurface at a large isocontour level with small volume tends to be significantly distorted on a coarse grid, so we still need to use a reasonably high grid resolution and a suitable C φ corresponding to a large enough volume for the assessment. For N = 256, we observe that the criterion R V > 0.5 is satisfied for R φ > 0.05, i.e., the corresponding data points are within the shaded region in Fig. 6.
As an illustrative example, we show the close-up view of VSF isosurfaces in the TG flow for Re = 400 at t TG = 7.0 with N = 512, N = 128, and N = 64 in Fig. 7. We observe two clear twisted vortex tubes from the reference numerical VSF solution with N = 512 and R φ = 0.5624 in Fig. 7a. By contrast, the vortex tube is barely resolved for the lowresolution VSF solution with N = 128 and R φ = 0.0508 in Fig. 7b. Although some parts are under-resolved, e.g., the very twisted part at the right end of the vortex tube in Fig. 7b, the essential tube structure is still preserved with R φ = 0.0508. On the other hand, the vortex tube is not well resolved in Fig. 7c with N = 64 and R φ = 0.0212. Therefore, we propose the criteria φ < 0.05 and R φ > 0.05 (17) to quantify both the VSF deviation and Lagrangian tracking performance. The threshold values are selected for the perturbed TG flow which has features of both turbulent and transitional flows with large-scale coherent and small-scale turbulent flow structures.

Morphology of vortex surfaces
We compare the morphology and dynamics of vortex surfaces in the flows under different magnitudes of the HIT disturbance. Figure 9 depicts the VSF isosurfaces with V φ = V /16 at two typical times t = 0.5 and 1.5 in cases M0, M1, and M2. The surfaces are color-coded by |ω|. The selection of V φ = V /16 is based on the satisfactory R φ in Fig. 8 and the visualization effect of displaying fine-scale structures. The transition from a smooth flow to a turbulent flow occurs around t TG = 4.5 ∼ 6.0 in TG flows, or equivalently t = 0 ∼ 1.5 in case M0. This process under the TG symmetries [21] has been elucidated by the evolutionary geometry and topology of VSF isosurfaces [6,11]. Initially, the two blob-like vortex surfaces approach to each other. They are flattened and then merged together around t TG = 4.5 or t = 0. As shown in Fig. 9a and b, the edge of the merged structure rolls up into vortex tubes during t = 0 ∼ 1.5. At late times, the tubes are tangled and twisted, constituting a turbulent-like flow. The VSF isosurface reveals the Lagrangian-like continuous evolution of vortex surfaces, which has essential differences from the visualization using Eulerian vortex criteria [5,6] such as Q [9] and λ 2 [10] criteria.
Next, we illustrate how the imposed HIT disturbance changes the evolution of vortex surface during the transition in case M1 with η M = 0.05 and in case M2 with η M = 0.1. First, the time evolution of φ and R φ in Fig. 8 indicates that the effective time period for the vortex surface tracking is around t = 0.5 ∼ 1.5 based on the criteria in Eq. (17). Compared with case M0, the imposed HIT disturbance in case M1 significantly wrinkles the vortex surface and produces more small-scale vortex tubes at t = 0.5 in Fig. 9b, and then some of the small-scale tubes are dissipated at t = 1.5 in Fig. 9e. On the other hand, the large-scale TG coherent structures are similar in cases M0 and M1. For more intense HIT disturbance in case M2, the vortex surface breaks down into a turbulent-like state with the emergence of numerous small-scale vortex tubes in Fig. 9c and f. The largescale TG structures are barely distinguished, while the regions of very large local vorticity magnitude remain due to the TG symmetries. Figure 10 provides a close-up view of the continuous evolution of vortex surfaces to reveal the local vortex dynamics, where the zoom-in region is marked by the black dashed   Fig. 9b in case M1 (η M = 0.05). The color bar is the same as that in Fig. 9 box in Fig. 9b. First, the perturbed vortex tube with large vorticity magnitude is wrinkled, then it is persistently stretched and twisted to form a spiral vortex tube, as marked in red boxes in Fig. 10a and f. By contrast, the wrinkled vortex surface with small vorticity magnitude tends to be smoothed as marked in blue boxes in Fig. 10a and f. The locally twisted spiral vortex tube reveals the important vortex core dynamics [31] in turbulent flows. This vortex geometry is consistent with the strained spiral vortex model [16] for small-scale turbulent structures, and supports the structure-based subgrid scale models [32] for the large-eddy simulation.
We continue to zoom in the region in the red box in Fig. 10 to show the detailed evolution of the twisted vortex tubes in Fig. 10. Besides the persistent twisting of vortex tubes with high |ω|, we notice that the vortex tubes are intertwined with each other in Fig. 11d,e,f. At the mean time, the twisted vortex tube is flattened into a ribbon-like structure, suggesting a possible way for the reduction of the characteristic length scale.
For comparison, Fig. 12 shows the VSF isosurfaces in the same regions in Figs. 10 and 11 in cases M0 and M2. The VSF isosurfaces in case M0 in Fig. 12a and b have much less wrinkling than those in case M1, and the surfaces in both the large and small |ω| regions stay smooth. The VSF isosurfaces in case M2 in Fig. 12c and d show much more wrinkled surfaces and chaotic small-scale structures with small |ω|. In particular, the comparison of Fig. 12b and d shows that the large-scale vortex tube structure in case M0 remains in case M2, and it consists of multiple small-scale, thread-like structures. This observation supports the strained spiral vortex model [16,33] in HIT. Moreover, the geometric differences of vortex surfaces in the three cases and the time coherence at   Fig. 9. The VSF isosurfaces with |ω| < |ω| max /10 in (a) and (c) and with |ω| < |ω| max /3 in (b) and (d) are cut off for clarity. Some vortex lines are integrated from points on the surface different times (e.g., t = 0.5 and 1.5) are not well characterized in the isosurface of |ω| (not shown).

Vortex surface distortion and local vorticity
From the visualization of VSF isosurfaces, the deformation of vortex surfaces tends to occur near the strong vorticity region. The local geometry of VSF isosurfaces is quantified by the mean curvature κ ≡ (∇ · n)/2 with n ≡ −∇φ/|∇φ|. As illustrated in Fig. 9, the surface distortion at the same time for cases M0, M1, and M2 increases with η M . The geometric difference between the TG flow and the perturbed TG flows is characterized by the volume-averaged mean curvature difference where κ M denotes the mean curvature in cases M1 or M2, and κ TG denotes the mean curvature in case M0. Thus, we quantify the preferential deformation in cases M1 and M2 using the volume averaged mean curvature difference conditioned on low and high vorticity regions 0 ≤ |ω|/|ω| max ≤ 1/8 and 7/8 ≤ |ω|/|ω| max ≤ 1 in case M0, respectively, in Fig. 13. It is clear that the surface distortion due to the HIT disturbance is enhanced within the high-|ω| region with time, while it is mitigated within the low-|ω| region, consistent with the observation in Fig. 10. This preferential surface distortion implies that the energy cascade from large to small scales tends to occur in high-|ω| region.

Conclusions
We propose criteria of tracking vortex surfaces in complex flows. The criteria in Eq. (17) involve not only the deviation φ of the numerical VSF solution, characterizing the attachment of vortex lines on VSF isosurfaces, but also the conservation ratio R φ , characterizing the Lagrangian-like tracking performance.
The new criterion on the numerical conservation is positively correlated to the volume overlap ratio between the isosurfaces of the numerical and reference VSF solutions. This criterion quantifies the error from the numerical artifact due to the nonuniqueness of VSF solutions in the correction step of VSF calculation. The comprehensive assessment of VSF solutions based on Eq. (17) provides a time range when the vortex surface tracking is satisfactory for both the accuracy of the vortex surface definition and the time coherence of the one-to-one flow map of VSF isosurfaces.
In order to assess the vortex tracking in chaotic but still coherent flow fields, we develop a turbulent-like flow combining the large-scale coherent structures in the TG flow and small-scale turbulent structures in HIT. The HIT disturbance is added to the TG flow near the TG transition time, which accelerates the transition to turbulence while preserves the large-scale TG symmetries. In the flow evolution, the vortex tracking criteria are satisfied within a finite period from t = 0.5 to 1.5 in cases M1 and M2 with HIT strengths η M = 0.05 and 0.1.
From the evolution of vortex surfaces during the effective tracking period, we demonstrate several local vortex dynamics. First, the imposed HIT disturbance significantly wrinkles vortex surfaces and accelerates the rolling up of small-scale vortex tubes. Second, the wrinkled vortex tube with large |ω| tends to be further twisted around the vortex Fig. 13 Temporal evolution of the averaged mean curvature difference in cases M1 and M2 conditioned on low-and high-|ω| regions axis under the self-induction, generating more small-scale structures and contributing to the energy cascade. By contrast, the wrinkling is mitigated in the region with small |ω|.
In future work, the tracking criteria can be applied to more complex flows, such as compressible turbulence [34], with more available computational resources. In particular, the detailed and quantitative investigation on the relation between the scale cascade of vortex surfaces and the energy cascade in spectral space is expected for high-Re HIT and wall-bounded turbulent flows.

Appendix A: Effect of the numerical scheme on the vortex surface tracking
The numerical regularization for Eqs. (6) and (7) is implemented using a dissipative finitedifference scheme for the convection term. For example, the WENO scheme is used in the present study. This scheme can effectively suppress numerical oscillations [26], while it may cause excessive numerical dissipation which reduces R φ and weakens the time coherence in the vortex surface tracking. Thus, we also test the eighth-order targeted essentially nonoscillatory (TENO) scheme [35,36], which is presumed to be less dissipative than WENO, for discretizing the convection term. We apply the WENO and TENO schemes to calculate the VSF in the TG flow for Re = 400 with different N and use the criteria in Eq. (17) to assess the performance of the two schemes for tracking vortex surfaces. The case parameters are listed in Table 2.
Compared with the WENO scheme, the TENO scheme has less numerical dissipation while has milder suppression on numerical oscillations. Figure 14 plots the temporal evolution of φ and R φ for the five TG cases. For the VSF deviation, both WENO and TENO perform well. In general, they have typically φ < 0.05 at t TG = 0 ∼ 6 with N ≥ 256, and φ for TENO is overall larger than that for WENO on the same grid. For the numerical conservation, TENO outperforms WENO, even when we compare the results for TENO on grid points N 3 and WENO on grid points (2N) 3 .
By comparing VSF visualizations from WENO and TENO schemes on the lowresolution grid with N = 128 in Fig. 15a and b, the TENO scheme clearly resolves the rolling up of vortex tubes, as the same as that from the reference solution in Fig. 15c. By contrast, the WENO scheme cannot resolve the vortex tube due to the excessive numerical dissipation in large straining regions. On the other hand, the TENO scheme can introduce more numerical oscillations in smooth regions and take more computational time than the WENO on the same grid (not shown). Hence, we still use the WENO scheme with a high enough grid resolution to calculate the VSF in the turbulent-like flows for the overall performance.  15 Comparison of the VSF isosurfaces in cases TG1, TG2, and TG5 with different numerical schemes and grid resolutions. The color bar is the same as that in Fig. 9