Vortex identiﬁcation methods applied to wind turbine tip vortices

. This study describes the impact of postprocessing methods on the calculated parameters of tip vortices of a wind turbine model when tested using Particle Image Velocimetry (PIV). Several vortex identiﬁcation methods and differentiation schemes are compared. The chosen methods are based on two components of the velocity ﬁeld and its derivatives. They are applied to each instantaneous velocity ﬁeld from the dataset and also to the calculated average velocity ﬁeld. The methodologies are compared through the vortex center location, vortex core radius and jittering zone. 5 Results show that the tip vortex center locations and radius have good comparability and can vary only a few grid spacings between methods. Conversely, the convection velocity and the jittering surface, deﬁned as the area where the instantaneous vortex centers are located, vary between identiﬁcation methods. Overall, the examined parameters depend signiﬁcantly on the postprocessing method and selected vortex identiﬁcation criteria. Therefore, this study proves that the selection of the most suitable postprocessing methods of PIV data is pivotal to 10 ensure robust results.


Introduction
The wake of a wind turbine is characterized by the presence of vortex structures. Two types of concentrated vortices can be identified, which are shed from the root and the tip region. The latter form strong helical shapes that influence the wake of the wind turbine. 15 The tip vortices are generated by the pressure difference between the top and lower side of the blade tip, which lead to a flow the pressure side to the suction side of the blade (Karakus et al., 2008;Sherry et al., 2013b). In this way, the tip vortices of a wind turbine represent a source of noise (Arakawa et al., 2005) and energy loss (Shen et al., 2005). Moreover, the wake development needs proper consideration in the layout of a wind park (Marten et al., 2020), as it can affect the performance of wind turbines located downstream. Therefore, a more detailed characterization of the wind turbine wake vortices does represent 20 a relevant research topic.
Since the first introduction of Particle Image Velocimetry (PIV) applied to wind turbine aerodynamics by Smith et al. (1990), a number of experimental investigations have been performed and a variety of methods have been employed to identify the vortex center and other characteristics. Yang et al. (2012) studied the formation and evolution of helical tip vortices of a wind turbine model under atmospheric boundary layer wind. It is shown, by using the vorticity to identify the vortices, a high variation in the position of the tip vortices. This effect is known as wandering or jittering and it is related to turbulence, vibrations of the model turbine (e.g. blades and tower) and the PIV system. Additional investigations (Maalouf et al., 2009;Soto-Valle et al., 2020) show the same effect using different identification methods such as the Q-criterion or circulation-based methods. Micallef et al. (2014) studied the mechanism of the initiation of the tip vorticity in a wind turbine. The findings showed how 30 the vorticity convects and forms a unique and symmetrical tip vortex behind the trailing edge. The location of the vortex center, identified by the maximum vorticity value, was found to be slightly inboard the rotor. Sherry et al. (2013a) studied tip vortices from a wind turbine in a water channel. The results highlighted the breakdown of the wake due to the mutual interaction between helical structures of the tip vortices, which is highly dependent on the tip speed ratio. Additionally, the jittering of the tip vortices was also detected. 35 Ostovan et al. (2019) studied the effect of winglets on the tip vortices. The results showed that the winglets cause faster downstream convection and a radially outward motion of the tip vortices, thus more wake expansion. The tip vortex was characterized by means of the velocity field when the convection velocity is subtracted, i.e. the vortex-induced velocity field, and the center was identified where the velocity is zero.
In summary, several vortex identification methods (VIM) have been employed so far. Yet, a common methodology for 40 the identification of vortices in the wake of a wind turbine has not been defined, as shown in Table 1. Furthermore, some investigations do not provide the complete implementation methodology, thus hindering an extensive comparison between methods. This paper aims at comparing different vortex identification methods to evaluate their suitability to study the tip vortices of 45 a wind turbine. The methods are applied to velocity field planes that were obtained through PIV in the near wake of a wind turbine model located in a wind tunnel facility (Soto-Valle et al., 2020). The comparison between methods and schemes is done by means of vortex center location, vortex core radius and jittering zone. In total three different VIMs are compared: vorticity, Q-criterion and Graftieaux. The first two VIMs require differentiation, thus, the application of six different schemes is examined.

50
The following section, Sect. 2, gives the mathematical overview of the methods to identify vortices and the differentiation schemes applied on their implementation. Subsequently, the wind tunnel and test rig used to generate the experimental dataset are introduced in Sect. 3, followed by the methodology in Sect. 4. The results are presented in Sect. 5 to conclude with the most important remarks in Sect. 6.

55
Many vortex identification methods have been proposed in the literature (Spalart, 1988;Hunt et al., 1988;Graftieaux et al., 2001;Vétel et al., 2010;Liu et al., 2016;Shkarayev and Kurnosov, 2017;Zhang et al., 2018a, b;Liu et al., 2020). In this work, three identification methods are compared. The chosen methods are based on the velocity field U , differing in their derivative orders. Consequently, the methods are based on the velocity field (U ) and first-order derivatives (∇U ). The selected methods are: 60 -Graftieaux's method, Graftieaux et al. (2001) -Vorticity magnitude, Spalart (1988) -Q-criterion, Hunt et al. (1988) Additional methods can be derived from the eigenvalue analysis, such as λ 2 , ∆ or swirling strength criteria (Zhang et al., 2018b). However, for the scope of this research, they represent similar approaches and therefore, only the selected methods are 65 analysed.
A full description of the selected methods is given below. For a more extensive review of VIMs, the interested reader is directed to Zhang et al. (2018b).

Graftieaux's method
This method identifies the vortex through a global quantity, Γ 1 , from an equivalent solid-body rotation. This function allows to 70 determine the location of the vortex center. Equation 1 shows the scalar, Γ 1 , defined in a discrete space.
where P is a fixed point to evaluate, U M is the velocity of the M surrounding points to P in the surface S, − − → P M is the radius vector that connects the point P with M , N is the total number of points considered in the surrounding of P , and z is the unit vector, normal to the surface plane S.
75 Figure 1 shows a graphic representation of the parameters for the calculation of Γ 1 . Over a two-dimensional frame, Γ 1 represents the topology of the surrounding flow to the point P . In this way, Γ 1 is the average contribution of the angles between the velocity U M and the radius vector. Therefore, at the vortex center, the value of Γ 1 tends to be close to the unity because the velocity contribution is perpendicular to the radius vector. Depending on the grid size, the value of Γ 1 might not reach unity as the center of the vortex could not be located on a grid 80 point. Therefore, the vortex center is estimated as the position of the maximum value of Γ 1 .

Vorticity magnitude
The vorticity is defined as the curl of the velocity field as show in Eq. 2 where ω is the vorticity. Equation 2 shows also a two-dimensional representation, where the out of plane vorticity is a function 85 of the velocity field U = U (u(x, y), v(x, y)) ≡ U (x, y).
In this way, the vorticity quantifies how the velocity vector changes when it moves in a direction perpendicular to it and therefore, is a natural candidate for vortex identification. Indeed, this method has been used for a long time (Spalart, 1988). In the vortex core, the vorticity is predominant compared to the shear rate deformation, due to the rotation of the fluid. Hence, the vortex core is identified as a region of high vorticity. At the vortex center, the vorticity reaches its maximum value. Con-90 sequently, the maximum value of vorticity can be used to locate the vortex center. However, as it has been pointed out by many authors (Liu et al., 2019;Zhang et al., 2018a), vorticity cannot distinguish between parallel shear motion and vortical motion. As an example, a laminar boundary layer shows high vorticity even though no vortical motion is present. Moreover, a vorticity threshold must be chosen to plot the vorticity iso-surfaces and the determination of the vortex core radius relies implement and commonly used in the literature.

Q-criterion
The most common vortex identification methods are based on the analysis of the velocity gradient ∇U . For instance, from the analysis of the eigenvalues of ∇U , Eq. 3, three invariants can be found (P, Q, R) In particular, the second invariant Q can be obtained through: where tr is the mathematical trace. In this way, the second invariant, Eq. 4, defines the Q-criterion. The method can be interpreted as the difference between the vorticity magnitude and the magnitude of the strain rate (Kolář, 2007). Hence, similar to the vorticity magnitude, the vortex core will be characterized by positive large magnitudes of Q since the rotation of the 105 fluid is predominant compared to the strain rate in this region. In addition, areas characterized by parallel shear motion, without rotation, will not be identified as a vortex (Q < 0) and therefore, overcoming one of the limitations of the use of vorticity in vortex identification. Nevertheless, a threshold is still needed in order to determine the core region.

Differentiation schemes
As shown in the previous section, several vortex identification methods are based on the gradient of the velocity field, then 110 inherently the evaluation of flow field derivatives is necessary. In this way, the differentiation of the velocity fields from either computational data or experimental techniques (as PIV) is needed. Both normally come in a discrete format.
In case of PIV data, the choice of the differentiation scheme becomes more relevant due to the presence of noise affecting the measurements. Noise sources include optical distortion, light sheet non-homogeneity, transfer function of CCD and particle characteristics, among others (Foucaut and Stanislas, 2002). Indeed, the process of differentiation can amplify the effects and 115 therefore compromise the results.
Many methods have been developed to calculate spatial derivatives from discrete data. The most frequently used methods are based on discrete differential operators applied to the surrounding points of the position to evaluate (Foucaut and Stanislas, 2002). In this way, the formulation of these schemes can be obtained through Taylor expansion. Equation 5 shows a generalization of the derivative scheme application on a function f over the dimension x in the point j (Raffel et al., 2018).
estimated from the uncertainty from the measurements (Lourenco and Krothapalli, 1995). Accordingly to Foucaut and Stanislas (2002), there is a trade-off between the truncation error and the noise amplification, therefore, the increment of the order will increase the uncertainty of the scheme.
The backward, central and forward differencing schemes provide the simplest implementation. Nevertheless, additional schemes have been studied with the purpose of either increasing the accuracy or reduce the uncertainty of the results. Raffel Another tested methodology shown in the same work is the least-squares scheme, a second-order scheme, designed to minimize noise propagation. However, this approach has the tendency to smooth the estimation of the derivative because the outer data is weighted more than the inner data.
The latter schemes are defined for a single variable function and applied in one dimension at a time. Conversely, the velocity 135 field can be influenced by the complete spatial coordinates. Therefore, the velocity gradient should depend on the surrounding flow. Raffel et al. (2018) proposed then, the circulation scheme that accounts for the effect of the surrounding flow (see Table   2). The first derivative is expressed as a central difference of derivatives in the other direction. This method reduces noise compared to the central difference scheme since the velocities of six neighboring points are considered instead of two. Table 2 shows a summary of the mentioned schemes with their accuracy and uncertainty.
140 Table 2. Summary of differentiation schemes and implementation a .

Experimental dataset
The analysis shown in the rest of this paper relies on the stereo-PIV-dataset from previous work by Soto-Valle et al. (2020). In the following, the experimental setup is presented.
The experiments were carried out in the closed-loop wind tunnel at the Technische Universität Berlin. The wind turbine 145 model, Berlin Research Turbine (BeRT) (Pechlivanoglou et al., 2015), was located in the large test section, with a crosssectional area of 4.2 × 4.2 m 2 . The freestream velocity was set to u ∞ = 6.5 ms −1 and the rotational speed to f = 3 Hz. Figure 2 shows the facility and test rig.
BeRT is a three-bladed, upwind horizontal axis wind turbine model with a rotor diameter of D = 3 m. The blades were twisted, tapered, and based on Clark Y airfoil profile along the full span. The Reynolds number was in the range of Re = 150 1.7 − 3.0 × 10 5 based on the chord length and the operational conditions.
The velocity field was measured using a stereo-PIV system, consisting of a Quantel Dual-Nd:Yag double laser with energy of 171mJ, a mirror arm, the laser sheet optics, two 14bit PCO 2000 cameras with a CCD-chip resolution of 2048×2048 pixel, and an ILA synchronizer connected to the measurement computer.
The measurement plane was horizontal and was centered at the tip location when the blade was in the horizontal position.

155
For the purpose of this study, one vortex age was chosen at φ = 40 • . This provides a total of 1200 recorded image pairs in a phase-locked position. The image postprocesing was done with the software PIVview3C (PIVTec GmbH), resulting in a field of view of 435 × 435 mm 2 and a spatial resolution of ∆x = ∆y = 3.6 mm. A two-dimensional analysis is carried out on the dataset. To conduct a three-dimensional analysis of the vortex structures, 160 additional parallel planes are needed. Therefore, only the two in-plane velocity components are used (x − y), even though the out of plane velocity w is available from the stereo-PIV measurements.
The application to obtain the vortex properties follows, while the statistical analysis of the available data is described at the end of this section.

Vortex center and convection velocity 165
The velocity fields are analyzed through the application of the VIMs described in Sect. 2. In the case of the vorticity magnitude and Q-criterion, VIMs are implemented using the differentiation schemes shown in  Results are presented in a normalized form and bounded by the unity. In the case of Γ 1 , the parameter already fulfill these requirements by its definition, while the results of vorticity and Q magnitudes are normalized according to Eq. 6: After the shedding from the tip, the tip vortex is both translating and rotating at the same time. Considering this, the convection velocity (downstream, x and outboard directions, y) is estimated as the velocity magnitude corresponding to the vortex center location. The latter is a common estimation in the literature (van der Wall and Richard, 2006;Yamauchi et al., 1999). 180 Therefore, the estimation is also affected by both the VIM and the scheme chosen on their application.

Core radius
The core radius is calculated using the following steps: 1. The induced velocity field U , Eq. 7, is obtained by means of the subtraction of the convection velocity from the velocity field. The resulting velocity field is characterized by presenting the induced contribution only.
2. The swirling velocity is analyzed through vertical and horizontal lines using the vortex center as an origin reference. The study is done in both directions to check the symmetry of the vortex, as vortices can have asymmetric shapes (Skinner et al., 2020). The procedure is repeated for each VIM and scheme and applied to both the average and the instantaneous velocity fields. Both vortex characteristics, convection velocity and core radius, are normally used to describe the evolution of the vortex at different ages (Snel et al., 2007;Nilsson et al., 2015) as well as a parameter to quantify induced drag penalties from the tip vortices (Ostovan et al., 2019).

Statistical analysis
A statistical analysis is made over the complete dataset of the instantaneous velocity fields. In this way the vortex center, the 200 convection velocity and the core radius are analyzed in terms of their location and magnitude variations.
An ellipse is used to define the jittering characteristic zone, similar to the work of Sherry et al. (2013b). The semi-axes of the ellipse a and b are defined to include all vortex center locations and the overall surface of the ellipse is calculated as S = πab.

Results
The results are presented as follows. First, an overview of the VIMs and schemes applied to the average velocity field are 205 shown using the vortex center locations and convection velocities based on the three identification parameters Γ 1 , ω, Q. The subscripts provide the information about the scheme implementation, e.g, ω CD shows the results from vorticity using the central difference scheme.
Second, a statistical analysis of the complete set of instantaneous velocity fields is performed. In this way, the location of the vortex center locations is studied in order to define the shape, distribution and surface of the jittering zone described by each 210 VIM and scheme. Additionally, these results are used to compare the scattering of the convection velocities and core radii. Figure 6 shows the magnitude distribution of the parameter Γ 1 after the application of Graftieaux's method with 8and 24points on the average velocity field. Both cases show a concentration of the magnitude in a core with one peak in the same location at (x/∆x, y/∆y) = (11, 14) and with an almost identical magnitude of Γ 1 ≈ 0.97, although the 24−p scheme extends 215 its distribution over a larger zone.

Average velocity field
Graftieaux's method has been developed for stationary vortices, while indeed, the test case is a superposition of the vortexinduced velocities and the streamwise flow, which convects the vortex downstream. The latter difficulty is overcome by subtracting the background velocity. Sherry et al. (2013b) proposed subtracting the average phase-locked velocity, U obtained by averaging each magnitude, streamwise (u) and lateral (v) from the full field of view.  On the other side, the presence of multiple maxima might also be due to small-scale structures within the vortex, as suggested by Bonnet (1998). It is conceivable that these structures are originated during the shedding of the tip vortex from the blade.
Certainly, the pressure difference between the pressure and suction sides of the blade is only one of many effects that take 2012). The involved structures are also affected by the blade shape, tip geometry, Reynolds number, and load distribution (Giuni and Green, 2013a) and generally merge into a single structure.
In conclusion, the multiple peaks could be caused by the uneven shedding of vorticity in the chordwise direction. In the work of Micallef et al. (2014) that shows the formation of the tip vortices in a horizontal axis wind turbine, a complex vorticity distribution along the blade chord is observed, which seems to cause multiple vorticity peaks inside the core. These multiple 245 peaks can be identified in the vortex core even after the tip vortex has been shed from the blade. In the present results, the same effect is obtained when the high uncertainty schemes are applied (CD, RE, BD and FD). Regarding the position of the vortex centers, the locations are shown in Fig. 9, overlapped with the vorticity magnitude distribution. It can be seen that the identification method does not have a strong influence on the estimation of the vortex center location with differences up to y/∆y = 2 and x/∆x = 4 grid steps in the lateral and streamwise directions, respectively. In case of vorticity and Q-criterion employing BD, CD, FD and RE schemes, the estimation of the vortex center is different from the geometrical center of the shape described by the core. As a result the convection velocity also differs significantly 255 between schemes. Figure 10 illustrates the axial and lateral velocities estimated from each VIM and scheme. Therefore, the estimation of the convection velocity is recommended with the smoother VIMs and schemes: Graftieaux or vorticity and Q-criterion while employing LS or CM schemes. Alternative methods such as the prediction from time series vortex locations might be also successful using the rest of the schemes due to the small discrepancy between the vortex center 265 locations between VIMs and schemes; however, more than one vortex age is needed.

Statistical analysis
The vortex center locations are identified on each instantaneous velocity field, which constitutes the complete PIV dataset. In the interest of clarity, only the Graftieaux 24-points as well only vorticity magnitude cases are presented. The reason is due to  points scheme applied to each velocity field. In fact, the zone can be highlighted as an ellipse that has its main axis on the lateral direction with 5 grid points more than the streamwise direction. Figure 11 b) shows the added up distribution along the streamwise direction. From Fig. 11, the jittering effect is clearly visible and agrees with previous results from fixed wings (Thompson, 1983;Giuni and Green, 2013b;Bandyopadhyay et al., 1991;Beresh et al., 2010), helicopter rotors (van der Wall and Richard, 2006; Mula et al., 2011) and wind turbines (Maalouf et al., 2009;Soto-Valle et al., 2020;Sherry et al., 2013a). According to these references, the source of the jittering can be varied such as geometry effects, wall boundary layer turbulence, free stream turbulence, surface irregularities or changes in the core structure. Additionally, the vibrations of either the model or the test 280 rig supports can produce small changes in the field of view, resulting in the meandering motion of the vortex. In this way, the jittering in Fig. 11 a) shows the spreading of the vortex center locations with y as a prevalent direction. The probability distribution over the streamwise direction, Fig. 11 b), fits very well with a normal distribution. These characteristics are in agreement with Sherry et al. (2013a), where it is shown that the jittering of the tip vortices in the wake of a horizontal axis wind turbine is predominant in the radial direction compared to the streamwise direction and that at early vortex ages, the 285 normal distribution is a good fit of the tip vortex center distribution. In agreement with the results presented here, Mula et al. (2011) also observed that the jittering of tip vortices generated by helicopter rotors present a preferential direction. Figure 12 shows the jittering zones for the Graftieaux method and vorticity calculations with the different differentiation schemes. For the purpose of clarity, only ellipse perimeters that contains the 100% of the vortex center locations are presented.
It is noticed that the ellipse described by the Graftieaux VIM is thinner in the streamwise direction observing approximately 290 two to four grid steps less than the vorticity VIM, depending on the schemes applied. In fact, the area described by the ellipses on the Graftieaux VIM is 80 ± 6 while in the case of the vorticity and Q are between 145 ± 24, corresponding to an 80% increment. The size variation between schemes in the vorticity VIM, Fig. 12, is more uniform in terms of directions and goes approximately to one grid step in any direction. Even though the area swept by the scattering of the estimated vortex center locations are similar in magnitude, in the case of 295 vorticity and Q, the probability distribution over these zones differs between the applied schemes. Figure 13 shows the contour diagram of two representative distributions applying vorticity VIM. Hence, Fig. 13 b), which is done with LS scheme, shows a more concentric distribution than CD scheme, Fig. 13 a). Q-criterion and the rest of the schemes exhibit similar results (see Least square scheme. Figure 14 shows the added up distribution towards the streamwise and lateral axes when the vorticity VIM and the CD 300 scheme are applied, together with fitting curves. Each fitting curve is chosen using higher the coefficients of correlation (R 2 ) between normal, binormal and Weibull distributions. In Fig. 14 a), it can be noticed that the spreading of the vortex centers is over two peaks with a distance of approximately 4 grid points and therefore, a binormal distribution fits better with the data with a coefficient of correlation of R 2 = 0.99. In case of the lateral direction, Fig. 14 b), a binormal and Weibull distributions exhibit good fitting, with the coefficients of correlations of R 2 = 0.97 and R 2 = 0.95, respectively.  Figure 15 shows the best-fit curves of the probability distribution when the differentiation schemes BD, CD and LS are applied. In the streamwise direction, Fig. 15 a), the LS scheme presents a normal distribution as the best fit. This is in agreement with the unique maximum observed in the analysis of the average velocity field. Instead, for all the other schemes the binormal distribution is the best fit. The difference is due to the smoothing properties of the LS scheme. In fact, the application is implemented using points up to the two positions farther in the grid, where the outer points are weighted more than the inner ones (see Table 2). Similarly, in the lateral direction, Fig. 15 b), the binormal distribution is likely the best fit with a coefficient of correlation between 0.97 − 0.99. In this case, peaks are very close, which also gives the Weibull distribution a good fit with coefficients of correlation between 0.95 − 0.97. In this direction, only one position is concentrated, consequently, the smoothing effect produced by the LS scheme does not have the same effect in this direction.

315
Overall, the estimation of the vortex center location is influenced by the VIM and scheme. In the same way, the convection velocity and the core radius are affected by the implemented methodology.
To see the effect of the VIM and schemes on the instantaneous convection velocity, Figure 16 shows the normalized convection velocity in the streamwise direction, u c /u ∞ . Each black dot represents the convection velocity estimated through the corresponding VIM and scheme on each velocity field from the PIV dataset. At the same time, the average of the dataset mag-320 nitudes is visible as a solid blue line. Moreover, the average value obtained when the methodology was applied on the average velocity field (see Fig. f:03 and Sect. 5) is displayed with a solid red line.
Graftieaux VIM, Fig. 16 a), presents a variation of 10% around its average value (blue line). In the case of the vorticity VIM, Figs. 16 b) and c), the variation increases up to 70% and 33% for the CD and LS schemes, respectively.
Several estimations fail, such as u c /u ∞ < 0 or u c /u ∞ > 1, as well in the average velocity analysis, due to the fact that they through the conditional averaging methodology. This error difference increases between 12 − 18% in the case of vorticity and Q-criterion (schemes CD, F D and RE) and to more than 50% with the scheme BD. The visible ring-like concentration on Fig. 17 contrasts the uneven shedding hypothesis from the average results and it leads to the idea of an artifact of the schemes, as it preserves the same ring-like structure even when an instantaneous velocity field is analyzed.
335 Table 3, shows the full set of average values and their corresponding standard deviations (SD). As discussed, vorticity and Q-criterion yield similar results. Nevertheless, the SD results from Q-criterion are equal or higher than vorticity in all the schemes and in both directions, which is presumable caused by the power of two on its derivatives implementation.
The average convection velocity is also estimated from the average flow field (red line). However, It should be noted that average results must not be overstated because of the tip vortex jittering (van der Wall and Richard, 2006). It is also remarkable 340 how convection velocity has the lowest results from BD scheme and on the contrary, the highest from FD scheme. In fact, both schemes ignore information either forward or backward from the grid on the implementation of differentiation. Therefore, they are not suitable for vortex analysis.  Figure 18 shows the normalized core radius, r c /∆x when it is estimated through along the horizontal axis. It is found that the results on the vertical axis are slightly higher, by 0.3∆x, which is small enough to assume the vortex as symmetric. In this way, the magnitude of the core radius is r c /∆x ≈ 3.4 independently on the VIM and scheme applied. The standard deviation varies between methods; however, it is always lower than 1∆x. This is expected due to the small number of grid 350 points that are found within the core (see Fig. 5). Indeed, the results from the average velocity field are also similar between VIMs and schemes but slightly higher than the statistical results, r c /∆x ≈ 4.4.
Several Vortex Identification Methods (VIMs) and implementation schemes have been applied to the two components of the velocity field data in the near wake of a wind turbine model, obtained through PIV measurements. The methodology was 355 applied to the average velocity field as well as the instantaneous velocities resulting in a statistical analysis of the PIV dataset.
In case of the average flow field, the chosen VIMs and schemes provide different magnitude distributions of the identification parameters. Nevertheless, the vorticity and Q-criterion yield the same estimations of the vortex center locations in all the schemes analyzed. Hence, as long as the vortex is well-formed the vorticity VIM is preferred over the Q-criterion because of the lower standard deviation results.

360
Through the statistical analysis, it is concluded that different methodologies lead to different interpretations of the tip vortex behavior. Even though the jittering zone is found to be ellipsoidal for all the VIM and schemes, the probability density function of the vortex center locations varies in the streamwise direction from one single peak with the Graftieaux, vorticity and Q (least-squares scheme) to a binormal distribution with the other implementations.
The two peaks found in the jittering are determined as an artifact produced by certain schemes. The latter can be avoided 365 using either Graftieaux VIM or vorticity and Q-criterion while employing the least-squares scheme.
It is concluded that the vortex center locations are within a small variation range and their comparability is viable independently on the VIM or scheme. Nevertheless, first-order schemes, such as backward and forward differences, should be avoided.
The convection velocity presented a higher dependency on the VIM and scheme applied. Therefore, and keeping in mind 370 that the results have shown good comparability regarding the vortex center locations, it is recommended to use the information of several vortex ages instead of the swirling velocity approach to estimate the convection velocity. Conversely, the vortex core radius only showed a grid step variation between VIM and schemes.
Overall, Graftieaux's method is the recommended VIM to track the tip vortex. Indeed, the method does not use differentiation and has shown to be independent of the number of neighboring points used. Moreover, it presents the lowest standard deviation 375 between all the methodologies applied here.
Data availability. Velocity field data and results can be provided by contacting the corresponding author.