Sensitivity Analysis of Maximum Circulation of Wake Vortex Encountered by En-Route Aircraft

Wake vortex encounters (WVE) can pose significant hazard for en-route aircraft. We studied the sensitivity of wake vortex (WV) circulation and decay to aircraft mass, altitude, velocity, density, time of catastrophic wake demise event, eddy dissipation rate, wing span, span-wise load factor, and WV core radius. Then, a tool was developed to compute circulations of WV generated/encountered by aircraft en-route, while disregarding unrealistic operational conditions. A comprehensive study is presented for most aircraft in the Base of Aircraft Data version 4.1 for different masses, altitudes, speeds, and separation values between generator and follower aircraft. The maximum WV circulation corresponds to A380-861 as generator: 864 and 840 m2/s at horizontal separation of 3 and 5 NM, respectively. In cruise environment, these WV may descend 1000 ft in 2.6 min and 2000 ft in 6.2 min, while retaining 74% and 49% of their initial strength, respectively. The maximum circulation of WV encountered by aircraft at horizontal separation of 3 NM from an A380-861 is 593, 726, and 745 m2/s, at FL200, FL300, and FL395, respectively. At 5 NM, the circulations decrease down to 578, 708, and 726 m2/s. Our results allow reducing WVE simulations only to critical scenarios, and thus perform more efficient test programs for computing aircraft upsets en-route.


Introduction
Wake vortices (WV) are produced by flying aircraft and may persist in the atmosphere for several minutes [1]. The hazards associated with WV encounters (WVE) have been thoroughly analyzed for aircraft on approach and departure phases [2][3][4][5][6][7][8], but for several reasons WVE en-route used to receive less attention, regardless the atmospheric conditions are often favorable for WV to remain strong for a long time, mainly because the natural atmospheric turbulence is generally low at cruise altitudes [9,10]. A first reason was that, when cruising, much altitude is available for recovery from strong WVE, which contributed to the perception that WVE do not pose a hazard. Secondly, the probability of severe WVE en-route is significantly lower than the probability of a similar situation when sequencing and merging arrival traffic flows in terminal airspace, or when giving take-off clearances in airports. Finally, since large WV-induced rolling deviations are extremely rare during cruise, and extreme vertical load variations can be attributed to many atmospheric irregularities such as clear air turbulence, a WV hazard was often believed not to exist at cruise [11].
Although the current rate of reported incidents related to WVE at en-route altitudes may be low, the problem of WVE en-route has received increasing attention in the last decade. This is due to the fact that WVE have become more frequent [1] and they will probably increase in the near future due to the expected evolution of the main factors contributing to WVE risk en-route (the characteristics of the generator and follower aircraft, encounter geometry, and tropopause altitude [12]) and other issues [13]: a growing amount of traffic in a very limited airspace (optimal cruise altitudes for the majority of jet aircraft lie in a rather thin layer around the tropopause); an increasing disparity in the size of

Modeling of Wake Vortex
The wake generated by aircraft forms a WV system composed of two counter-rotating vortices. According to the Kutta-Joukowsky theorem, for a generator aircraft in level flight, their initial circulation Γ 0 is: where g is the gravity acceleration, ρ is the air density, and m, U ∞ , and b are the mass, TAS, and wing span of the generator aircraft, respectively, and s is the span-wise load factor: the ratio between the initial lateral spacing between the vortices b 0 and b. The vortices sink due to the mutually induced velocity, and also experience lateral motion (called transport), depending on the prevailing wind. The lift distribution and mass of the generator aircraft have a significant influence on this trajectory. Particularly, the initial WV is formed about 10 wing spans behind the aircraft and then starts descending and decaying [21], i.e., the WV circulation decreases with time, starting from Γ 0 . The decay and descent rates depend mainly on the wind shear, atmospheric turbulence, and thermal stratification [10,22,23]. The latter is the atmospheric condition that causes the greatest effect on the WV evolution, and is frequently described by the Brunt-Väisälä frequency N [12,24]: where θ is the potential temperature, ω is the potential temperature lapse rate, and h is the altitude. Higher N is associated with higher stability and buoyancy force acting on the WV, which causes faster WV decay rate [12]. N is generally higher in the stratosphere than the troposphere. Moreover, an increase in the tropopause altitude due to global warming has been reported [25]. This is concerning since the probability of WVE is usually larger in the troposphere [1], i.e., from sea level (SL) to the tropopause, usually at around 11 km at medium latitudes. In addition, the severity of the WVE increases for increasing tropopause altitude [12]. Thus, works such as the present research, focused on WVE in en-route phase, will become increasingly important. Unfortunately, there is a lack of validated WV decay models for cruise flight. This is why Luckner and Reinke [21] used Sarpkaya's decay model although it has only been tested against measurements at low altitudes, leaving uncertainty regarding its applicability at cruise altitudes. Namely, researchers claimed that additional efforts should be directed at improving WV models for en-route applications after observing that the WV descent rate was in one case probably much higher than predicted by the model [12]. Anyway, the decay in Sarpkaya's model affects only the effective WV circulation Γ [26]: where t is time and t c is the time of the catastrophic wake demise event. This model does not account for the effect of stratification on WV decay, which can be neglected for low stratification levels of the atmosphere [21], as usual below the tropopause [10]. Conversely, the eddy dissipation rate (EDR) ε has a significant effect on t c , given a normalized EDR ε * and a normalized time of catastrophic wake demise event t c * : Aerospace 2021, 8, 194 4 of 22 the relation between both variables is t c * = 9.18 − 180ε * , for ε * < 0.0121 or which is most interesting for our work, given the typical values of ε * [27]: For modeling the decay, researchers have also used a Betz WV model based on far field conservation principles, modified with an empirically based core size [28], and the probabilistic/deterministic two-phase WV transport and decay model (P2P/D2P), describing WV decay and transport based on atmospheric conditions [29][30][31]. In particular, P2P/D2P describes the vortex decay and descent through two consecutive decay phases, following large eddy simulation (LES) results [32,33]. In the first phase (the diffusion phase), the normalized circulation Γ * = Γ/Γ 0 as a function of the non-dimensional time t * is [31]: while in the second phase (the rapid decay phase), it is [31]: In the P2P/D2P model, Γ is the average over circles of radii from 5 to 15 m, Γ 5−15m . The characteristic decay parameters (i.e., A, R * , T * 1 , ν * 1 , T * 2 , and ν * 2 ) are shown in [31] as a function of atmospheric conditions (e.g., N and ε), where applicable (see details later on). The dependency of T * 2 on ε * for null N * for the mentioned LES data is also shown, together with a model by Sarpkaya relating ε * with T * 2 for null N * . The velocities induced by the WV on the flow field can be obtained by superimposing the two single vortices (left and right vortices, which are assumed to have identical circulation but of opposite direction [21]), using the tangential velocity V t model of Burnham-Hallock [34], based on [35], which yields good results [21,[36][37][38]: where r is the radial distance from the WV center line and r c is the WV core radius. The WV sink rate w WV can be computed based on Equation (10) by imposing r = b 0 [21]: Integrating Equation (11), we obtain the WV altitude descent (or sinking) h WV with time [21]:

Sensitivity Analyses
Sensitivity analyses were done on the effect on WV circulation, decay and sinking of variations of several parameters [20]. From the propagation theory, the sensitivity of a parameter y to variations of a parameter x is: If we analyze the sensitivity of Γ 0 to changes in generator aircraft mass m, it appears from Equations (1) and (13) that a variation of a given order of magnitude in m causes a Aerospace 2021, 8, 194 5 of 22 variation of Γ 0 of the same order of magnitude. However, if we consider that, in levelled rectilinear horizontal flight, the aircraft lift L is: where S is the wing lay-out area and C L is the lift coefficient, then, changes in m would propagate into changes of Γ 0 half as large: A previous work reports that variations in m propagate not directly into exactly the same variations of Γ 0 but into slightly smaller variations [21]. We arrive to similar results if in the previous equations we assume a linear relationship between S and m as in [39], i.e., S = a I + a I I m = 19.463 + 1.645m, obtained using data from [40] (a similar approach is followed in [36], where a functional relation is established between the MTOW and various aircraft parameters, such as S, based on data of existing aircraft): The parameter A in Equation (17) is 0.5 for m = 0 and tends to 1 with increasing value of m (see Figure 1); e.g., A > 0.90 for m > 47.5 tons, following the linear correlation between S and m [39]. The key is that Equation (16) would be the one relevant to ATCO and ATM (interested in the sensitivity of Γ 0 to the mass of an aircraft already built, i.e., with a given value of S). Conversely, Equation (17) would be relevant to aircraft manufacturers for aircraft design.
A previous work reports that variations in propagate not directly into exactly the same variations of Γ 0 but into slightly smaller variations [21]. We arrive to similar results if in the previous equations we assume a linear relationship between and as in [39], i.e., = + = 19.463 + 1.645 , obtained using data from [40] (a similar approach is followed in [36], where a functional relation is established between the MTOW and various aircraft parameters, such as , based on data of existing aircraft): The parameter in Equation (17) is 0.5 for = 0 and tends to 1 with increasing value of (see Figure 1); e.g., > 0.90 for > 47.5 tons, following the linear correlation between and [39]. The key is that Equation (16) would be the one relevant to ATCO and ATM (interested in the sensitivity of Γ 0 to the mass of an aircraft already built, i.e., with a given value of ). Conversely, Equation (17) would be relevant to aircraft manufacturers for aircraft design. as in [39].
As per the sensitivity of Γ 0 to changes in flight altitude ℎ, from Equations (1) and (18), we can obtain Equation (19): where is the adiabatic coefficient for air and ′ is the universal gas constant divided by the air molecular mass. If we consider that and follow the International Standard Atmosphere (ISA) [41], while all the other parameters are constant with altitude since the error is negligible (e.g., the change in from SL to ℎ = 20 km is below 0.63%), then: Figure 1. Parameter A vs. generator aircraft mass m, assuming linearity between wing lay-out area S and m as in [39].
As per the sensitivity of Γ 0 to changes in flight altitude h, from Equations (1) and (18), we can obtain Equation (19): where γ is the adiabatic coefficient for air and R is the universal gas constant divided by the air molecular mass. If we consider that ρ and T follow the International Standard Atmosphere (ISA) [41], while all the other parameters are constant with altitude since the error is negligible (e.g., the change in g from SL to h = 20 km is below 0.63%), then: where T SL and ρ SL are the ISA temperature and density at SL (288.15 K and 1.225 kg/m 3 , respectively), and k is the ISA temperature lapse rate (−6.5 K/km). In the lowest layer of the stratosphere, from 11 to 20 km, k = 0 K/km, thus T = T 11 is constant, and: where T 11 and ρ 11 are the ISA temperature and density at the tropopause (h 11 = 11 km). Thus, the sensitivity of Γ 0 to the flight altitude h can be finally expressed as: where and are the ISA temperature and density at SL (288.15 K and 1.225 kg/m 3 , respectively), and is the ISA temperature lapse rate (-6.5 K/km). In the lowest layer of the stratosphere, from 11 to 20 km, = 0 K/km, thus = 11 is constant, and: where 11 and 11 are the ISA temperature and density at the tropopause (ℎ 11 = 11 km). Thus, the sensitivity of Γ 0 to the flight altitude ℎ can be finally expressed as: From Equation (23), for example, an altitude increment of 1000 ft at flight level FL195, i.e., ∆ℎ/ℎ = 5.13%, causes an increase of Γ 0 of 3.78%, while an increment of 1000 and 2000 ft at FL460, i.e., ∆ℎ/ℎ = 2.17% and 4.35%, causes an increase of Γ 0 of 4.80% and 9.60% (see in Figure 2). If we consider the sensitivity of Γ 0 to ∞ , a variation of a given order of magnitude in ∞ causes a variation of Γ 0 of the same order of magnitude but with opposite sign (the same occurs if we consider variations of density , Mach number , span-wise load factor , wing span , or initial lateral spacing between vortices 0 ): However, in levelled rectilinear horizontal flight, from Equations (1) and (14), we obtain Γ 0 = ∞ /(2 0 ), hence: and the same occurs if we consider variations of wing lay-out area or lift coefficient . Regarding the decay in effective WV circulation Γ, from Equations (3) and (13), a variation If we consider the sensitivity of Γ 0 to U ∞ , a variation of a given order of magnitude in U ∞ causes a variation of Γ 0 of the same order of magnitude but with opposite sign (the same occurs if we consider variations of density ρ, Mach number M, span-wise load factor s, wing span b, or initial lateral spacing between vortices b 0 ): However, in levelled rectilinear horizontal flight, from Equations (1) and (14), we obtain Γ 0 = U ∞ SC L /(2b 0 ), hence: and the same occurs if we consider variations of wing lay-out area S or lift coefficient C L . Regarding the decay in effective WV circulation Γ, from Equations (3) and (13), a variation of a given order of magnitude in Γ 0 causes a variation of Γ of the same order of magnitude. Moreover, based on the previous sensitivity analyses for Γ 0 , we can deduce the sensitivity of Γ to the parameters that affect Γ 0 . As per the sensitivity of Γ to t c , from Equation (3) we obtain this expression (see example plot of C in Figure 3):  Regarding the sensitivity of Γ to changes in , since affects the value of , as shown in Equations (4) to (7): For 0.0121 < * < 0.2535, using Equation (6) we obtain (see example plot of in Figure 4 left): • For 0.2535 < * , using Equation (7)   Regarding the sensitivity of Γ to changes in ε, since ε affects the value of t c , as shown in Equations (4)-(7): • For 0.0121 < ε * < 0.2535, using Equation (6) we obtain (see example plot of D in Figure 4 left): • For 0.2535 < ε * , using Equation (7) we obtain (see example plot of D in Figure 4 right): Aerospace 2021, 8,194 8 of 22 • For 0.2535 < * , using Equation (7) (27) and (right) Equation (28).
When analyzing the WV sink rate w WV , from Equations (3) and (11): a variation of a given order of magnitude in Γ 0 causes a variation of w WV of the same order of magnitude. Based on the previous sensitivity analyses, we can deduce directly the sensitivity of w WV to t c and the parameters that affect Γ 0 , with the exceptions of b, r c , s, and b 0 , which are now treated joining Equations (1) and (29): Taking advantage of the typical relationship between the WV core radius r c and b (i.e., r c is usually defined as a small percentage of b [36][37][38][42][43][44]), the sensitivities of w WV to variations in b, r c , s, and b 0 are as follows: Thus, changes in b and r c propagate into changes of w WV twice as large in absolute terms. Almost the same occurs for s and b 0 : using typical values of b, r c (i.e., 1% to 5% of b [36][37][38][42][43][44]), s (i.e., 0.75 to 0.85 [3,36]), and b 0 , E yields values above 1.990. Particularly, the higher the values of b 0 , b, or s, or the lower the value of r c , the closer is E to 2 (see Figure 5). However, as long as r c sb (which is mostly the case), a change of r c does not impact the WV descent. What matters most for this research is the hazard of large aircraft, where the mutual velocity induction is compared to the region of potential vortex (e.g., in P2P/D2P, b 15 m and we use Γ 5−15m ). Thus, the impact that matters is the variation of b 0 due to wing design (inboard and/or outboard loading) and dynamic variations terms. Almost the same occurs for and 0 : using typical values of , (i.e., 1% to 5% of [36][37][38][42][43][44]), (i.e., 0.75 to 0.85 [3,36]), and 0 , yields values above 1.990. Particularly, the higher the values of 0 , , or , or the lower the value of , the closer is to 2 (see Figure 5). However, as long as ≪ (which is mostly the case), a change of does not impact the WV descent. What matters most for this research is the hazard of large aircraft, where the mutual velocity induction is compared to the region of potential vortex (e.g., in P2P/D2P, ≫ 15 m and we use Γ 5−15 ). Thus, the impact that matters is the variation of 0 due to wing design (inboard and/or outboard loading) and dynamic variations of 0 en-route for new aircraft (e.g., the A350 uses the flaps for having variable camber during flight).

Methodology
In our research, a circulation generator module (CGM) was developed, which produces circulation values and other outputs that can be used as input for computing the upsets in follower aircraft. As an example, the following values of initial WV circulation were used in previous works: 343, 683, and 783 m 2 /s [21], and 530 m 2 /s [24]. However, the CGM is able to compute the maximum possible circulation of WV generated or encountered by aircraft for a given set of conditions, while disregarding unrealistic scenarios. Among other benefits, this allows keeping to a minimum the costly computations of upsets in follower aircraft to search for the most severe upsets. The main input data necessary to execute the CGM are:

1.
Generator aircraft and/or follower aircraft type: Any aircraft for which the necessary data are available in BADA version 4.1 [17] can be used; 2.
Flight Level (FL) of the corresponding aircraft: In a previous work on WVE, most of the tests were made at FL370 (37,000 ft), and the influence of WVE altitude was studied at FL410 [21]. In our case, the CGM is able to find the maximum circulation either for a specific FL or in a given range of FL, from a minimum specified FL (e.g., FL200, as this is a typical lower limit for cruising FL) to the highest possible FL (i.e., the ceiling of the corresponding aircraft for the given flight conditions); 3.
Generator-follower separation: This must be provided in terms of horizontal distance or time (e.g., 5 NM, as in [21]), to compute the decay of the encountered WV and/or vertical separation.
Once these input data are introduced, the CGM computes the initial circulation Γ 0 of the WV downwind of the generator aircraft using Equation (1). This can be done for a range of FL, masses, and speeds of the generator aircraft (see details below), generating multiple scenarios with their corresponding Γ 0 . Then, the CGM disregards unrealistic scenarios, i.e., impossible combinations of mg, U ∞ , and FL for the generator and follower aircraft, based on their performance. For this purpose, aircraft performance data from BADA is used. For example, an aircraft may generate WV with very large Γ 0 for a given FL, but the follower aircraft may not be able to reach that FL, for any of the scanned combinations of mg and U ∞ . In brief, the method is based on computing the maximum rate of climb (ROC) that the aircraft can achieve for the given flight conditions. If, for those conditions, that value is below a minimum specified ROC (see details below), the corresponding scenario is discarded. Then, among the non-discarded scenarios, the maximum Γ 0 is obtained. Finally, Γ is computed using Sarpkaya's and the P2P/D2P decay models, and the maximum circulation of the WV encountered by the follower aircraft is established based on the input aircraft separation. In all the computations done in this work, the following hypotheses were assumed: 1.
Atmospheric properties: The values of the atmospheric variables (e.g., ρ, temperature T, pressure p, kinematic viscosity ν, and the speed of sound a) depend on the FL, as established by the ISA [41]; 2.
Eddy dissipation rate (EDR): EDR ranging from 10 -8 to 2 × 10 -7 m 2 /s 3 and 10 -6 to 10 -1 m 2 /s 3 appear in [24,45], while Meischner et al. [46] measured a maximum of 0.05 m 2 /s 3 inside strong storm cells. The CGM allows using an EDR value determined by the flight altitude, as in [47] or setting a given EDR value at the user's discretion. In our case, the EDR was set to 10 -6 m 2 /s 3 , assuming neutral stratification (i.e., null N) and low levels of atmospheric turbulence, as in [21,45,48]. This is natural at just above/below the tropopause [10,24], and is a conservative approach, since those atmospheric conditions imply the largest hazard potential of WVE [12], i.e., this low EDR corresponds to worst-case scenarios as per potential upsets on the follower aircraft, since the WV decay is slow and the WV persistence is high; 3.
Aircraft mass: In a research on WVE, aircraft weights in steps of 5% ranging from 65% to 95% of the maximum take-off weight (MTOW) were considered [12]. In another work with the A330-300 aircraft, three masses were used [21] WV core radius or size: Usually, the initial WV core radius r c0 is specified between 1% and 5% of b [42][43][44]. A smaller r c for a given Γ 0 is a conservative approach for hazard considerations [49], as confirmed by our sensitivity analysis around Equation (32): smaller r c leads to faster sink rate. From flight test measurements, r c0 is 3.5% of b [36][37][38]. While r c grows with increasing vortex age, it does not change much (e.g., the decay in Sarpkaya's model does not affect r c ), and it does not change either with the atmospheric conditions. Hence, a constant WV core radius r c = r c0 was assumed in [36]. Likewise, we kept r c constant at 3.5% of b (4.5% of b 0 ) from BADA. Different simulations revealed that r c has no significant effect on the upsets on follower aircraft [50,51], while a previous work [8] suggests otherwise; 6.
WV span-wise load factor: This factor s depends on combined effects of the loadings of the wing and horizontal tail plane [3], and usually ranges from 0.75 to 0.85. For the wing, it is often assumed an elliptical lift distribution [3], which is particularly valid for aircraft flying en-route above FL285 [12]. For these wings, s is commonly assumed to be π/4 [36]. Again, s does not change much with time or the atmospheric conditions (except maybe with stratification), particularly, assuming low levels of turbulence, as long as the decay rates are moderate, s is constant [21]. Thus, we assumed that s is constant at π/4; 7.
Minimum rate of climb: In most of the controlled airspace, a minimum ROC is imposed by ATM regulations and practices. Thus, optimal cruise altitudes are typically the highest FL that allow for this minimum ROC, given the aircraft mass and maximum thrust that can be achieved for that particular altitude (note, however, that wind conditions aloft might significantly change this cruise altitude, leading in some cases to lower optimal altitudes with more "benign" winds). This minimum ROC is typically 500 ft/min [52], and this value is used in this work to disregard unrealistic combinations of m, U ∞ , and FL. Nevertheless, this parameter is at the user's discretion, and might be configured to any other value; 8. Turbofan thrust model: BADA includes performance data for various turbofan thrust control models. For computing the maximum ROC, we used the maximum climb (MCMB) turbofan rating model; 9.

Practical Application: Maximum Circulation of Wake Vortex Generated by Aircraft
In a first analysis, the CGM was configured to compute the maximum possible circulation of WV generated by multiple aircraft in realistic scenarios. The input data were: Generator-follower separation: We used the current minimum horizontal separation en-route in radar environment (i.e., 5 NM), and vertical separations of 10 and 20 FL (1000 and 2000 ft). Seeking for worst-case scenarios to test the feasibility of reduced separation standards, horizontal separations of 0.5 and 3 NM were also tested (the rationale would be the possibility of increasing airspace capacity, based on the requirement of standard deviation of lateral path-keeping error of 0.3 NM [55]). Table 1 shows the obtained results for several generator aircraft, and their flight operation conditions. Here, Γ 0 is the initial WV circulation, and Γ Sarp and Γ D2P are the circulations of the WV from Sarpkaya's and the P2P/D2P decay models, respectively, at different horizontal separation distances d SEP between the generator and follower aircraft. The table also shows the time separation t SEP between the aircraft for the given flight velocity U ∞ , and the WV sinking h WV , equivalent to the generator-follower vertical separation. Namely, the maximum values of Γ Sarp obtained in this analysis, at horizontal separation of 0.5, 3, and 5 NM, were 896, 864, and 840 m 2 /s, respectively, corresponding to the A380-861 as generator. The equivalent results for Γ D2P were 859, 829, and 805 m 2 /s (see also Figure 6). Table S1, provided as supplementary material, shows the results obtained for all the aircraft for which the necessary data are available in BADA version 4.1.  [26] and P2P/D2P [31] decay models, and data for * = 0.01 and * = 0.00, taken from [24], for A380-861 aircraft.
Several simple tests allowed checking the robustness of the CGM and the consistency of the results. Firstly, aircraft weights were sampled up to MTOW, on one side, and up to 95% of MTOW, on the other. As expected, this did not have an effect on the obtained maximum WV circulations, since scenarios with MTOW at cruise altitudes are discarded Figure 6. Wake vortex circulation Γ vs. time t from Sarpkaya's [26] and P2P/D2P [31] decay models, and data for ε * = 0.01 and N * = 0.00, taken from [24], for A380-861 aircraft. Table 1. Maximum circulations of the wake vortex generated by several aircraft and flight conditions of the generator. Several simple tests allowed checking the robustness of the CGM and the consistency of the results. Firstly, aircraft weights were sampled up to MTOW, on one side, and up to 95% of MTOW, on the other. As expected, this did not have an effect on the obtained maximum WV circulations, since scenarios with MTOW at cruise altitudes are discarded by the tool as unrealistic. Secondly, Mach values were sampled up to 97% of the M mo , on one side, and up to the M mo , on the other. Again, as expected, this did not have an effect on the obtained maximum WV circulations, meaning that these are always obtained not at the highest but the lowest tested Mach numbers (i.e., the M mrc in this case). This is explained by Equation (1) and the fact that, for given conditions in levelled rectilinear horizontal flight, flying at lower U ∞ requires using higher values of C L , as shown in Equation (14), and thus higher intensity of the wing tip vortices. Despite this, we did not test Mach numbers lower than M mrc since they are not interesting for nominal cruise operations, as the aircraft would fly slower and with higher fuel consumption (and potential future linear holding operations, where these low Mach values may be interesting [56], are out of the scope of this work).

Generator Mass (t) U ∞ (m/s) M (-) FL (-) r c (m) b 0 (m) Γ 0 (m 2 /s) d SEP (NM) Γ Sarp (m 2 /s) Γ D2P (m 2 /s) t SEP (s) h WV (t SEP ) (ft)
Finally, as shown in Table 1 and Figure 7, the overall maximum circulation is generally obtained at the ceiling of the corresponding aircraft (the ceilings were obtained using a criterion described later on, and thus may not be coincident with public aircraft specifications or the maximum operating altitude in BADA). The relation between flight altitude and circulation in levelled rectilinear horizontal flight is indeed interesting: Equation (14) alone would not provide a conclusive answer, since, to keep the lift constant, the decrease of air density with altitude could be compensated by either increasing U ∞ or C L . For given flight conditions, increasing C L would lead to higher intensity of the wing tip vortices (i.e., higher Γ 0 ). On the other hand, Equations (1) and (14) combined show that: Aerospace 2021, 8,194 13 of 22 m). However, the WV strength by P2P/D2P in the second decay phase is significantly lower than that given by Sarpkaya's model, and reaches 0 much earlier. Finally, the results from Misaka et al. [24] are higher than those from the P2P/D2P and Sarpkaya's models in the first half of the first decay phase, and lower in the second half. For instance, the difference between Misaka's results and P2P/D2P is on average −13.8%, with a standard deviation of 24.8%. For the same case mentioned above, Figure 8 shows the WV descent in altitude vs. time, as obtained from [24], the Burnham-Hallock model and Sarpkaya's circulation de- Thus, if we increase U ∞ or the altitude, or we decrease density, Γ 0 would increase (this agrees with the assessments made around Equations (23)- (25), which provide highly valuable insights). However, there are a few exceptions for which the overall maximum circulation is achieved at intermediate FL, not the ceiling or close to it. This is likely because the aircraft mass also affects Γ 0 , as shown in Equation (1). That is, these exceptions would be cases in which the drop in Γ 0 due to the lower flight altitude is counterbalanced by the fact that the aircraft can fly at those lower FL with much higher mass than at the ceiling.
For the case in which the A380-861 generates the maximum WV circulation, Figure 6 shows the circulation decay, as well as data from [24]. The circulation at generator-follower aircraft horizontal separation of 3 and 5 NM, and 1000 and 2000 ft, i.e., 304.8 and 609.6 m, below the altitude of the generator aircraft is also shown. P2P/D2P usually leads to vortices with longer lifetimes than Sarpkaya's model in stratified environments [29], but shorter lifetimes for neutral stratification (N ≈ 0, as we assumed). This is why the WV lifetime by Sarpkaya's model is longer in all our tests, as shown in Figure 6. P2P/D2P may predict slightly higher (up to 6.2% higher) circulations than Sarpkaya's model around the transition from first to second decay phase, but in general the results from both are very similar in the first phase (difference of 1.3% on average with a standard deviation of 3.1%, the −4.2% gap of P2P/D2P with respect to Sarpkaya's model in the first instant of time is due to the fact that P2P/D2P uses a circulation averaged over circles of radii from 5 to 15 m). However, the WV strength by P2P/D2P in the second decay phase is significantly lower than that given by Sarpkaya's model, and reaches 0 much earlier. Finally, the results from Misaka et al. [24] are higher than those from the P2P/D2P and Sarpkaya's models in the first half of the first decay phase, and lower in the second half. For instance, the difference between Misaka's results and P2P/D2P is on average −13.8%, with a standard deviation of 24.8%.
For the same case mentioned above, Figure 8 shows the WV descent in altitude vs. time, as obtained from [24], the Burnham-Hallock model and Sarpkaya's circulation decay, shown in Figure 6. The WV sinking at generator-follower aircraft horizontal separation of 3 and 5 NM (50 and 83 m, in 22 and 37 s, respectively), and a descent of 1000 and 2000 ft below the generator aircraft altitude are also indicated. Note that, in cruise level environment, the WV of an A380-861 may descend 1000 ft in 2.6 min and 2000 ft in 6.2 min, while retaining a significant fraction of its initial strength: 74.4% and 48.8%, respectively, according to Sarpkaya's model or 72.4% and 36.8%, according to P2P/D2P. In addition, it is worth recalling that Sarpkaya's model tends to predict faster weakening of the WV descent rate than P2P/D2P, resulting in smaller overall WV descent [29]. This agrees with Misaka et al. [24], who claim that WV may descent more than 2500 ft, well beyond the 1000-ft vertical aircraft separation prescribed in RVSM airspace, and typical en-route WV persistence of 2-3 min and sink rate of 400 ft/min [1]. In one case, researchers observed a WV descent rate probably much higher than predicted by the WV sink rate model [12]. Hence, encounters with WV of even higher circulations could happen well below the altitude of the generator aircraft. . In addition, it is worth recalling that Sarpkaya's model tends to predict faster weakening of the WV descent rate than P2P/D2P, resulting in smaller overall WV descent [29]. This agrees with Misaka et al. [24], who claim that WV may descent more than 2500 ft, well beyond the 1000-ft vertical aircraft separation prescribed in RVSM airspace, and typical en-route WV persistence of 2-3 min and sink rate of 400 ft/min [1]. In one case, researchers observed a WV descent rate probably much higher than predicted by the WV sink rate model [12]. Hence, encounters with WV of even higher circulations could happen well below the altitude of the generator aircraft. Figure 8. Wake vortex descent in altitude ℎ vs. time from the Burnham-Hallock model [34], and data for * = 0.01 and * = 0.00 [24], for the A380-861 aircraft.

Practical Application: Maximum Circulation of Wake Vortex Encountered by Aircraft
In a second analysis, the CGM was configured to compute the maximum possible circulation of WV encountered by several aircraft in realistic scenarios. The input data were: 1. Generator and follower aircraft type: F100-650, A320-212, A330-301, B772LR, and A380-861 were used in this analysis (a total of 25 combinations of generator-follower); 2. FL of the follower aircraft: Three different FL were considered: FL200, the ceiling of the aircraft for the given input conditions, and an intermediate FL; 3. Generator-follower separation: We used the current minimum horizontal separation en-route in radar environment (i.e., 5 NM), and vertical separations of 10 and 20 FL (1000 and 2000 ft). Seeking for worst-case scenarios to test the feasibility of reduced separation standards, horizontal separations of 0.5 and 3 NM were also tested (the rationale would be the possibility of increasing airspace capacity [55]). Table 2 shows, for the scanned masses and Mach numbers for the A320-212, realistic combinations of flight operation conditions at FL200, FL300, FL395, and FL398 (its ceiling). Tables A1, A3, and A5 in the Appendix A show, respectively, for the A330-301, B772LR, and A380-861, realistic combinations of flight operation conditions at FL200, FL300, and their ceilings. Table 3 shows, for several generator-follower horizontal and vertical separations, the maximum circulations of WV encountered by any aircraft that can fly at the indicated FL, behind a selected set of generator aircraft. The parameters shown are the same as in Table 1. The overall maximum circulation obtained in this analysis from  [34], and data for ε * = 0.01 and N * = 0.00 [24], for the A380-861 aircraft.

Practical Application: Maximum Circulation of Wake Vortex Encountered by Aircraft
In a second analysis, the CGM was configured to compute the maximum possible circulation of WV encountered by several aircraft in realistic scenarios. The input data were:
FL of the follower aircraft: Three different FL were considered: FL200, the ceiling of the aircraft for the given input conditions, and an intermediate FL;

3.
Generator-follower separation: We used the current minimum horizontal separation en-route in radar environment (i.e., 5 NM), and vertical separations of 10 and 20 FL (1000 and 2000 ft). Seeking for worst-case scenarios to test the feasibility of reduced separation standards, horizontal separations of 0.5 and 3 NM were also tested (the rationale would be the possibility of increasing airspace capacity [55]). Table 2 shows, for the scanned masses and Mach numbers for the A320-212, realistic combinations of flight operation conditions at FL200, FL300, FL395, and FL398 (its ceiling). Tables A1, A3 and A5 in the Appendix A show, respectively, for the A330-301, B772LR, and A380-861, realistic combinations of flight operation conditions at FL200, FL300, and their ceilings. Table 3 shows, for several generator-follower horizontal and vertical separations, the maximum circulations of WV encountered by any aircraft that can fly at the indicated FL, behind a selected set of generator aircraft. The parameters shown are the same as in Table 1. The overall maximum circulation obtained in this analysis from Sarpkaya's model, at horizontal separation of 3 and 5 NM, was 593 and 578 m 2 /s at FL200, 726 and 708 m 2 /s at FL300, and 745 and 726 m 2 /s at FL395, corresponding to the A380-861 as generator aircraft (the results at FL398 are not shown as they are very close to FL395). As expected, for the reasons explained before, the obtained maximum WV circulations increase with altitude. Table 2. Realistic flight operation conditions for the A320-212 at FL200, FL300, FL395, and FL398.  Table 3. Maximum circulations of wake vortex encountered by aircraft behind several generator aircraft in realistic scenarios, and flight conditions of the generator.  Tables A2, A4 and A6 in the Appendix A show, respectively, the ceilings of the A330-301, B772LR, and A380-861, the same parameters as in Table 3. Note that, if the generator and follower aircraft are flying at the same FL with separation of 3 or 5 NM, the WV sinking is around 100 or 200 ft, respectively. For 5 NM, the sinking is large enough so that the follower may likely fly above the WV, considering the typical error in altitude. Conversely, if the separation is reduced to 0.5 NM, the follower may perfectly encounter the WV left behind by a generator at its same FL (altitude-keeping errors from systems such as flight data recorders show that height-keeping errors below 50 ft predominate, while the proportion of errors beyond 300 ft is less than 2.0 × 10 −3 [55]).

Conclusions
The sensitivity of wake vortex (WV) to several parameters was studied: (a) changes in aircraft mass m propagate into changes of initial circulation Γ 0 half as large or almost as large, depending on m and whether the wing lay-out area S is given or we assume linearity between S and m; (b) the higher the flight altitude, the more sensitive is Γ 0 to altitude variations; (c) a variation of a given order of magnitude in true airspeed U ∞ , density ρ, Mach number M, span-wise load factor s, wing span b or initial lateral spacing between vortices b 0 causes a variation of Γ 0 of the same order of magnitude; (d) a variation of a given order of magnitude in Γ 0 causes a variation of Γ and WV sink rate w WV of the same order of magnitude; (e) finally, changes in b and r c propagate into changes of w WV twice as large in absolute terms, and almost the same occurs for s and b 0 . Then, a tool named circulation generator module (CGM) was developed to compute the circulation of WV generated/encountered by aircraft in en-route scenarios. The tool disregards impossible combinations of aircraft mass, true airspeed and altitude, based on aircraft performance data from BADA version 4.1. The maximum possible WV circulations were computed with the CGM in a comprehensive set of en-route scenarios, scanning several parameters for multiple generator and follower aircraft. The main conclusions are: Since experimental data are not available yet for validation purposes, flight tests could greatly help assess the accuracy of the proposed approach and the obtained results; and 5.
The CGM and our results can be used to assess the severity of WVE in the en-route phase. Indeed, the worst-case WV circulations obtained under realistic operational conditions allow reducing WVE simulations only to critical scenarios, and thus to perform more efficient test programs for computing aircraft upsets.    Data Availability Statement: Some or all data, models or code used during the study were provided by a third party (Base of Aircraft Data (BADA) version 4.1). Direct requests for this material may be made to the provider as indicated in the Acknowledgments. Some or all data, models or code that support the findings of this study are available from the corresponding author upon reasonable request (the circulation generator module (CGM)).

Acknowledgments:
The authors would like to thank, on one side, the valuable feedback from Ramon Dalmau about the information in Base of Aircraft Data (BADA), and, on the other side, EUROCONTROL, for giving us the license to use the BADA version 4.1. Direct requests for this material may be made to EUROCONTROL.

Conflicts of Interest:
The authors declare no conflict of interest.

Notation
The following symbols are used in this paper:

Appendix A
Further results from the CGM are included in this appendix as well as in Table S1 as supplementary material. Particularly, Tables A1, A3 and A5 show, respectively, for the A330-301, B772LR, and A380-861, realistic combinations of flight operation conditions at FL200, FL300, and their respective ceilings (FL415, FL427, and FL431), for the scanned masses and Mach numbers. Tables A2, A4 and A6 show, for several generator-follower horizontal and vertical separations, the maximum circulations of WV encountered by any aircraft that can fly at the ceilings of the A330-301, B772LR, and A380-861, respectively. The parameters are the same as in Table 3. Table A1. Realistic flight operation conditions for A330-301 at FL200, FL300, and FL415.