Far field wake vortex evolution of two aircraft formation flight and implications on young contrails

Abstract Large-eddy simulations (LES) have been employed to investigate the far-field four-vortex wake vortex evolution over 10min behind an aircraft formation. In formation flight scenarios, the wake vortex behaviour was found to be much more complex, chaotic and also diverse than in the classical single aircraft case, depending very sensitively on the formation geometry, i.e. the lateral and vertical offset of the two involved aircraft. Even though the case-by-case variability of the wake vortex behaviour across the various formation flight scenarios is large, the final plume dimensions after vortex dissolution are in general substantially different from those of single aircraft scenarios. The plumes are around 170 to 250m deep and 400m to 680m broad, whereas a single A350/B777 aircraft would produce a 480m deep and 330m broad plume. Formation flight plumes are thus not as deep, yet they are broader, as the vortices do not only propagate vertically but also in span-wise direction. Two different LES models have been employed independently and show consistent results suggesting the robustness of the findings. Notably, 
$CO_{2}$
 emissions are only one contribution to the aviation climate impact among several others like contrails and emission of water vapour and nitrogen oxides, which would be all affected by the implementation of formation flight. Thus, we also highlight the differences in ice microphysical and geometrical properties of young formation flight contrails relative to the classical single aircraft case.

Unterstrasser (1) UG2014 Unterstrasser and Görsch (2) Greek Symbols x, y, z mesh size in span-wise, flight and vertical direction, m t numerical time step, s eddy dissipation rate, m 2 /s 3 t=0 initial magnitude of circulation, m 2 /s 1 , 2 , 3 , 4 circulation of the vortices from left to right, m 2 /s σ init width of Gaussian plumes of FAC, m θ 0 MGLET potential temperature ω y vorticity perpendicular to y-axis, s ω vorticity magnitude, s

INTRODUCTION
Migratory birds flying in a flock improve their aerodynamic efficiency, save energy and increase their range (3,4) . Empirical evidence is given by heart rate records of white pelicans where Weimerskirch et al. (5) find that individuals flying at the back of a flock have lower heart rates and hence reduce their energy expenditure. Similarly, formation flight (FF) can increase the performance in the civil and military aviation sector. In close FF, aircraft (AC) have stream-wise separations of a few wing spans. Due to safety issues, however, close FF is only relevant for the military sector. For the commercial sector, extended formation flight (with separations of 10 to 40 wing spans) is a viable option as the risks of collisions become acceptable (6,7) . Follower AC benefit from flying in the up-wash region created outboard of a leading AC. Hereby, the induced drag is reduced, the lift-to-drag ratio increases and fuel consumption is lower as demonstrated in numerous numerical, wind tunnel and real flight studies e.g. (8)(9)(10)(11)(12)(13)(14) . These studies treat formations of two or three AC and aim at determining the sweet spot, i.e. the relative position of follower aircraft for which the lift-to-drag ratio is highest or the engine thrust is lowest. Roughly speaking, fuel benefits of around 10% can be expected for such formations.
The drag reduction is spatially inhomogeneous across the follower AC's wings and induces a rolling and pitching moment. Kless et al. (15) shows that trimming the rolling and pitching moments by aileron deflections reduces the drag benefits by around 10% (of the aforementioned 10% saving) in transonic conditions. Apart from such deflections of the control surfaces, differential engine thrust between left and right engine or an asymmetric fuel load in the wings can be alternative trim mechanism. Employing a combination of the two latter trim mechanisms, Okolo et al. (16) found that the full benefit of formation flight can be fully retained.
So far, aerodynamic benefits and resulting fuel savings during an actual formation flight were regarded. However, re-routing of and coordination among several AC is required to establish a formation for certain segments of their flight routes. Clearly, the fuel savings during the formation flight must substantially outweigh the re-routing induced fuel penalty. Considering a large network of solo mission flight routes, it is a complex optimisation problem to find candidates which should optimally join and build formations. The larger the network of cooperating airlines and the number of possible routes to be joined, the higher is the potential reward. Xu et al. (17) find net fuel burn reductions of nearly 8% and 6% for a 150-flight Star Alliance transatlantic schedule and a 31-flight single airline long-distance schedule, respectively. Corresponding reductions in direct operating costs are smaller because re-routing involves longer travel times.
The wake vortices of the lead AC may be advected by cross-winds. In particular in extended FF, their lateral displacement may not be negligible. In such cases it is not sufficient to simply keep the relative positions of the aircraft, but to account for the real position of the wake vortices. Finding the sweet spot, then, involves the sensing and detection of the wake vortices by the follower AC. In the best case, the flight path of the follower AC is autonomously adjusted by sophisticated flight control systems. Various real-time wake tracking algorithms have been proposed and tested recently (18,19) . They rely on wing-distributed on-board pressure sensor data.
Despite these challenges formation flight could be introduced with less technological efforts and adaptations compared to other efficiency options like natural laminar wings, blended wing bodies or open rotor techniques. Hence, this topic should receive attention by the scientific community as well as by aviation stakeholders.
In the context of greener aviation, fuel benefits translate into reduced CO 2 emissions. Notably, CO 2 emissions are only one contribution to the aviation climate impact among several others like contrails and emission of water vapour and nitrogen oxides (20,21) . The contrail radiative forcing (RF) is probably larger than the RF of the total accumulated CO 2 emissions from aviation (22,23) . Contrails consist of ice crystals that grow in moist environments e.g. (24,25) . They may live for many hours (26,27) and undergo a transition into contrail-cirrus, i.e. contrails which lose their line shape and resemble naturally formed cirrus e.g. (28) . If several contrails are produced in close proximity, they compete for the available atmospheric water vapour and mutually inhibit their growth (29) . Such a contrail cluster is expected to have a smaller climate impact than the total effect of the same number of individual contrails that grow isolated from other contrails. The properties of contrail-cirrus depend on atmospheric conditions, but also on initial processes in the aircraft wake e.g. (2,25,30) . The early contrail evolution is dominated by the wake vortex descent and ice crystal loss due to adiabatic heating in the descending primary wake e.g. (31,32) (the primary wake is the part of the wake that moves downward with the wake vortices) and is traditionally simulated by LES models coupled to an ice microphysical model (1,2,33,34) .
The present study deals with the evolution of the wake vortex system in formation flight scenarios. Large-eddy simulations (LES) are initialised with a four-vortex system of a two-AC formation and cover its evolution and decay over 10min. To our knowledge, such simulations have not been described elsewhere in the literature. The paper will highlight that the wake vortex evolution is much more complex than in the classical case behind a single aircraft. In order to check the robustness and to increase the fidelity of the LES results treating such an unprecedented case, two different LES codes, EULAG-LCM and MGLET, were employed separately. Moreover, the study details the implications on the properties of young contrails. A follow-up numerical study will focus on the contrail-cirrus evolution over several hours and assess the climate benefits by contrail saturation effects of formation flight scenarios. The present paper is divided into four sections. Section 2 describes the LES models and their setups. Section 3 presents the EULAG-LCM simulations results on the wake vortex and contrail evolution. The implications are discussed and summarised in Section 4. The appendix completes the work with numerical sensitivity tests and a comparison with MGLET wake vortex simulations.

NUMERICAL MODELS AND SETUP
Two LES models are separately employed in this study. Both codes have a long tradition of wake vortex simulations of single aircraft.
Employing EULAG-LCM, Unterstrasser et al. (35) investigates the dispersion of a passive aircraft exhaust tracer at cruise altitudes under the influence of the downward moving wake vortex pair. Unlike to MGLET, the finite-difference dynamical solver EULAG is additionally equipped with the Lagrangian ice microphysics code LCM which allows to simulate the contrail life cycle and the interplay of ice microphysics and wake vortex dynamics (1,2) .
MGLET is a finite-volume code, and its temporal version was successfully applied to study the dynamics of wake vortices after roll-up until decay. With this approach, a vortex pair with a constant velocity profile along flight direction is initialised. It may incorporate different vortex profiles and even multiple vortex systems. This technique enables taking into account various atmospheric conditions like turbulence, thermal stability and wind shear, both out of ground (36) , as well as in ground proximity (37,38) . Wake vortex simulations with MGLET were compared to water towing tank experiments (39) . Here, especially the interaction of the vortices with ground is well captured by the simulations. In Stephan et al. (40) , lift and drag coefficients retrieved from numerical simulations with MGLET were compared to flight tests. Additionally, qualitative agreement with lidar measurements at Vienna Airport was demonstrated. Main advantages of MGLET over EULAG are its fourth-order compact spatial scheme treatment and the inclusion of a dynamic sub-grid scale model.
The results section is entirely based on EULAG-LCM simulations. MGLET simulations are presented in the appendix and serve as a benchmark reference in the model comparison.
Throughout the paper, the transverse (cross-stream) direction is along x, the flight (streamwise or axial) direction along y, and the vertical direction along z, as illustrated in Fig. 1.
The following subsection shortly describes EULAG-LCM, whereas a MGLET-description is deferred to the appendix.

LES model EULAG with Lagrangian ice microphysics LCM
The base model EULAG (41,42) solves the Navier-Stokes (NS) momentum and energy equations. In its current version, various approximations of the NS equations can be dealt with, ranging from a Boussinesq, an anelastic, a pseudo-incompressible to a compressible version (43) . It relies on the iterative upwind scheme MPDATA (Multidimensional Positive Definite Advection Transport Algorithm (44) ) with substantially reduced implicit diffusion compared to the classical upwind scheme. The transport algorithm belongs to the class of NFT (non-oscillatory forward-in-time scheme) and is at least of second order. EULAG is an allscale model with atmospheric applications ranging from global (45) , mesoscale (46) to local (47) scale optionally featuring dynamic grid deformation (48) and adaptive moving meshes (49) . Subgrid turbulence modelling is based on a classical TKE approach (50,51) . Domain decomposition in all three spatial directions allows perfect scalability in massive parallel applications with O(10 4 ) processors (52) . In the present study, we use the anelastic formulation, a time-constant uniform Cartesian grid, and 2D horizontal domain decomposition with 32 × 6 processors. The EULAG solution procedure, the underlying equation set and numerical parameter choices used here are summarised in Section 2.1 of Unterstrasser et al. (35) . Fully coupled to the dynamical solver is the Lagrangian ice microphysics model LCM (53,54) , which comprises explicit aerosol and ice microphysics for simulating pure ice clouds like natural cirrus (55) or contrails (56,57) . In the Lagrangian approach, ice crystals are represented by simulation particles (SIPs) that are advected by the fluid. Each SIP represents a certain number of (identical) real crystals and stores information, e.g. on the discrete position x SIP , the mass m SIP and the habit of the ice crystals that are bundled in a SIP. In its full version, microphysicsrelated processes on the simulation particles include homogeneous freezing of liquid supercooled aerosol particles, heterogeneous ice nucleation, deposition growth of ice crystals, sedimentation, aggregation, latent heat release and radiative impact on particle growth. For the given problem, it is not necessary to use the full LCM apparatus. As in previous EULAG-LCM studies of the early contrail evolution (1,2) , deposition growth (including a Kelvin correction term) and latent heat release are the only microphysical processes switched on.
LCM uses information of the velocity and thermodynamic variables as provided by EULAG. A second order Runge-Kutta scheme is used to solve the advection equation for each SIPẋ where the particle velocity u SIP is a superposition of the fluid velocity u LES , an auto-correlated turbulent contribution u SGS , which is based on the TKE-value provided by EULAG and accounts for SGS motions, and lastly the terminal settling velocity w sed in the vertical direction (as written above, sedimentation is not relevant in the present problem and w sed is set to zero).
Microphysical processes are computed for each SIP, e.g. the ice crystal deposition growth is a function of temperature, water vapour and ice crystal propertieṡ Summing upṁ SIP of all SIPs in a grid box, one can derive the rate of change of the water vapour concentration qv, which appears as source term in the EULAG prognostic equation of qv LES . Similarly, latent heating as derived in LCM is accounted for in the EULAG temperature equation.
From the SIP data, Eulerian representations of, e.g. ice crystal number concentration N can be derived a-posteriori. Moreover, those 3D fields can be integrated along one or two spatial dimensions in order to obtain transverse/vertical profiles of ice crystal number (N t and N v ) or the total ice crystal number per meter of flight path (N ). Note that the derived quantities are averages along flight direction. Following U2014, we use the following definitions: EULAG-LCM simulations were run at the HPC cluster of DKRZ Hamburg. Typically, the massively parallel computations use 192 cores and run for 25h consuming around 5,000 CPUh.

Basic parameters
We perform temporal LES assuming homogeneous flow conditions along flight direction. In our previous wake vortex studies, the initial flow field consisted of a background turbulent field (produced by an a-priori simulation with specified turbulence level and stratification) and a pair of counter-rotating Lamb-Oseen vortices. Here, a four-vortex system is superimposed on the background field and describes the flow state downstream of the follower aircraft (FAC). This initialisation plane can be seen in the left part of Fig. 1. The vortices of the FAC are labelled V 1 and V 2 , and the ones of the leader aircraft (LAC) V 3 and V 4 . We will refer to V 2 and V 3 as inner vortices (IV) and to V 1 and V 4 as outer vortices (OV). The idealised 2D flow field of the initialisation plane is uniformly initialised at each slice of the simulation domain as indicated in the right part of the figure. We note that the black boxes are only used for illustration purposes and the actual 2D flow field is prescribed over a larger cross-section than suggested by the black boxes.
For the sake of simplicity, we study a formation flight scenario with two identical aircraft. We use two aircraft of type A350 or B777 with a wing span of b s = 60.9m and a mass which would generate wake vortices with circulation t=0 = 520 m 2 /s. The separation distance between the vortices is b 0 = (π/4) b s = 47.3m.
In a first simplistic and clearly unrealistic approach, each vortex V i is initialised with | i | = t=0 . A vortex core radius of r c = 4m is prescribed. We will refine the flow field initialisation in Section 2.2.2. Figure 1 sketches the geometry of the formation flight scenario. The two AC have a stream-wise separation DY (along flight direction) and span-wise (lateral) separation DX. In commercial applications the extended formation flight with stream-wise separations of N y = 20 to 40 wing spans is the preferred pattern and corresponds to a time separation of 5-10s (simply given by t LAC/FAC = N y b s /U, where U is the AC speed). Due to safety reasons, shorter separations are foreseen only in the military sector. Upon the passage of the FAC, the vortex pair of the LAC descends. In the plane of initialisation, which is downstream of FAC and illustrated by the shaded area in Fig. 1, the separation DY along flight direction is translated into a vertical displacement DZ of the LAC vortex pair (DZ = w 0 t LAC/FAC , where w 0 = (2 0 )/(π 2 b s ) is the descent speed of the vortex pair). A lateral (span-wise) separation DX of around 0.8 wing spans was found to give optimal fuel benefits in previous studies and is used in this study as a default. The red blobs indicate the relative position of the four vortex centres (VC) and the arrows show the sense of rotation.
The temperature at cruise altitude is T * = 217K, the relative humidity with respect to ice is everywhere RH * i = 110%, the Brunt-Väisälä frequency is N BV = 1.15 · 10 −2 /s and the eddy dissipation rate is = 10 −7 m 2 /s 3 . The pressure at the bottom of the domain is 250hPa. EULAG simulations use the anelastic approximation and a stable background temperature profile is prescribed according to N BV . MGLET uses the Boussinesq approximation and we prescribe θ 0 = 330K and d θ s /dz = 4.46K/km. The ambient conditions represent typical conditions of the upper troposphere and have been chosen in analogy to Unterstrasser (1) (abbr. from now on as U2014). Note that unlike to the lower atmosphere, where cloud formation is initiated at water vapour saturation, the formation of ice clouds in the upper troposphere is inhibited and substantial supersaturation is required for ice nucleation to proceed. This makes ice supersaturation and RH i -values > 100% a common phenomenon of the upper troposphere. Contrail ice formation occurs in the expanding exhaust jets and is not bound to ambient RH i . By the way, this gives also the explanation for skies that are crowded by contrails in an otherwise cloud-free scenario.
Finally, we summarise contrail-related aspects of the initialisation, which are similar to the A350/B777-setup of Unterstrasser and Görsch (2) (abbr. from now on as UG2014). The green discs in Fig. 1 show the positions of the exhaust plumes. At the time of initialisation, the plumes of the LAC are assumed to be fully entrained into the wake vortices and discs with uniform concentrations are collocated with the VCs (default setup of UG2014). For the FAC plumes, Gaussian plumes are initialised inboard of the VCs (setup as in Section 3.4 of UG2014, a more detailed description of this setup type is given in U2014). The initial ice mass and crystal number per meter of flight path are I 0 = 3.0 · 10 −2 kg/m and N 0 = 6.8 · 10 12 /m (total of both AC). This corresponds to an 'emission' index EI iceno = 2.8 · 10 14 per kg of burned fuel. Each plume carries one fourth of all ice crystals. Overall, around 82 · 10 6 Lagrangian particles are used to represent the ice crystals. As in preceding studies the initial ice crystal sizes are log-normally distributed with width r SD = 3.0 (details, see Section 2.2 of U2014).
The model domain has dimensions L x = 768m, L y = 132m and L z = 600m. The cruise altitude of the two AC is at z = 350m. Periodic boundary conditions are applied in the horizontal and rigid boundaries in the vertical. In several plots the z-coordinate will be shifted to identify the cruise altitude with z = 0m. Moreover, the x-coordinate will be shifted such that the body Numerical Parameters x, y, z 1m, 2m, 1m Ice crystal parameters of the (left) FAC is aligned with x = 0. Figure 1 illustrates the domain size (yet for a grid sensitivity simulation with L y = 792m) and the coordinate system. All basic parameters are summarised in Table 1.

Inner vortices
So far, we assumed that the two IVs are not disturbed and 2 = − 3 = t=0 . The fuel saving of the FAC results from the fact that it receives lift from the IV of the LAC. The lift loading is not equally distributed across the two wings of the FAC and | 2 | < | 1 | (15) . Moreover, the roll-up process of V 2 is affected by the close-by LAC IV V 3 . On the other hand, also the evolution of the LAC IV V 3 is disturbed as the vortex or certain areas of it impinge on the FAC wing. Such interactions between a solid body and a stream-wise oriented vortex have been simulated recently. A series of studies analysed the flow field a few chords downstream of a rigid flat plate wing with a specified angle-of-attack and results were found to depend sensitively on the lateral and vertical offset of the incident vortex relative to the wing (58,59) . Barnes et al. (60) considered a flexible wing including aeroelastic effects. The experimental study of Inasawa et al. (61) addresses the vortex impingement on a trailing wing by visualising the flow field with PIV (particle image velocimetry). They found that the circulation of the emerging wing tip vortex of the trailing wing can be higher or lower compared to the default case without an impinging vortex, depending on the relative position of the impinging vortex. All On the other hand, the findings may be generalised to extended FF scenarios as the properties of the lead vortex do not change too much between a few wing spans (once the roll-up is completed) and 30 wing spans behind the aircraft. However, those studies do not provide information whether or not the flow of the LAC IV V 3 is also disturbed far from the centre and whether or not the roll-up of the FAC IV V 2 is such that a potential vortex is a good approximation for the flow field at large radii. All in all, this leads to the situation that the flow field initialisation has some highly uncertain aspects.
In our default setup, we leave the two OVs V 1 and V 4 as is (using a standard Lamb-Oseen radial profile of tangential velocity u θ ) and damp the two IVs. For radii r < R weak , u θ is multiplied by a scaling factor that increases linearly from 0.2 at r = 0 to 1 at r = R weak . The modified profile u θ (r) together with the standard profile is plotted in Fig. 2(a). The consequences on the radial distributions of vorticity ω y (r) and circulation (r) are illustrated in panels (b) and (c).
Based on observations by McKenna et al. (62) we introduce additional turbulent velocity fluctuations (white noise fluctuations with root mean square velocity u rms,ET = 2 m/s) in a rectangular box encompassing the two IVs (see yellow box in Fig. 1). The box is L x,ET = DX broad and the lateral boundaries are aligned with the centres of the AC bodies. The vertical boundaries are 15m above the centre of V 2 and 15m below the centre of V 3 giving a box height of L z,ET = DZ + 30m.
Keeping in mind the uncertainties of the IV initialisation, several options are tested in an extended sensitivity experiment. In one simulation, both IVs are damped over a larger area by increasing R weak from 15 to 25m; see the red dashed lines in Fig. 2 for the impact on the vortex properties. The overly simplistic approach outlined in Section 2.2.1 can be referred to as simulation with R weak = 0m. The centres of the two IVs are around 15m apart from each other. Hence, a variation of R weak substantially affects the interference among the IVs.
In further tests, the vortex strength is reduced not only locally around the core region, but globally for all radii. This means that the circulation values | 2 | and | 3 | are decreased. In such cases with 'full' damping, the damping is uniform across the whole radius range and no additional local damping is introduced (i.e. R weak = 0m). In one simulation, both IVs are similarly weakened by setting 2 = − 3 = 0.5 t=0 . In another simulation, an unsymmetrical damping with 2 = 0.6 t=0 and 3 = −0.4 t=0 is prescribed. This assumes that the  destruction of or energy drain from the LAC vortex V 3 is stronger than the potential beneficial reduction of the induced drag at the FAC, which is linked to generation of a weaker vortex V 2 . In the full damping cases, not only the mutual interference is affected, also the effect of the IVs on the two OVs is modified. Finally, in another sensitivity simulation with default R weak = 15m, local turbulence is not enhanced. All sensitivity simulations concerned with the IV initialisation are listed in blocks 4 to 6 of the 'Sensitivity experiments' Section of Table 2. The sensitivity to the IV initialisation will be reported in Section 6.1.

Sensitivity experiments
Besides sensitivity tests accounting for the uncertainties of the IV initialisation discussed in the latter section, several physical parameters are varied. In real flight, it will not be operationally feasible to keep the lateral and vertical offset fixed over time, in particular the transverse shift has to consider the effect of cross winds. To account for these uncertainties, DX takes values of 45, 55 and 60m in a simulation series, complementing the default value of 50m. Moreover, DZ is reduced to 0 or 7m or enlarged to 25m (default 14m). A dominant parameter for the contrail evolution is the ambient relative humidity RH i , which is varied from 100% to 140% similar to previous studies. A list of all sensitivity simulations is given in Table 2. We note that the default values are framed by a box and that in each sensitivity simulation only a single parameter is varied.

RESULTS
The results section starts with a short presentation of an example simulation. This is followed by a detailed analysis of how the four vortices interact with each other and how they move. Finally, the implications on plume dimensions and contrail properties are discussed and differences to the classical single aircraft case are highlighted.

Example simulation
In this paragraph we shortly describe the flow/vortex evolution of an example simulation over several minutes. For this, Fig. 3 displays 3D contour plots of vorticity magnitude (the time and the contour surface level are indicated in each panel). We choose the default simulation, yet for illustration purposes with an enlarged domain length along flight direction (L y = 792m instead of 132m, the simulation is listed in the last row of Table 2). In the first panel at t = 2s, the four vortex tubes are straight and the relative positions between them reflect the geometric initialisation as outlined before in Fig. 1. Moreover, coloured contours are plotted in four slices perpendicular to the flight direction. Those reveal the enhanced turbulence that was superimposed in a rectangular box around the two IVs during the initialisation. The panel for t = 16s shows that secondary vorticity structures (SVSs) have emerged around the two IVs. Moreover, those two vortices move towards the right OV. Section 3.2 will explain in detail the reasons for this lateral propagation. Soon after, the right OV V 4 (as labelled in Fig. 1) captures V 3 . Those two now form a vortex pair and move away from the other IV V 2 . After 60s, the vortex pair V 3 &V 4 effectively moved upwards, many SVSs shape up around them and high vorticity values are not any longer confined to two concentrated tubes (as is obvious from the contour slice in front). On the other hand, the vortices V 1 and V 2 still feature straight vortex tubes. A closer inspection shows that vortex V 1 is stronger than V 2 , as the red area with vorticity magnitude ω > 1.5/s is larger. Moreover, a small-amplitude small-wavelength meandering of V 2 is apparent. After 3min, the vortex pair V 3 &V 4 has dissolved and only a few patches with elevated ω-values occur. Contrarily, the two other vortices can be still tracked as they still feature fairly straight vortex tubes and only weak SVSs are present. After 5min, the weaker V 2 -vortex has dissolved and only V 1 remains. Another 1.5min later at t = 6.5 min, mostly only irregular patterns of enhanced vorticity have survived. Most simulations are carried out until 10min, where we can safely assume that wake vortices induced dynamics has ceased.

Tracking of the vortex centres
In the classical case of a single aircraft, a pair of counter-rotating vortices rolls up and the vortices mutually induce a downward movement at an initial speed of At cruise conditions, the final vertical displacement and the decay of the vortices is mainly governed by the strength of thermal stratification e.g. (35,63) and the well-known Crowinstability may occur (64) . The formation flight simulations show a much more diverse behaviour and their early evolution is strongly affected by an interplay of the two close-by IVs. This will be demonstrated next by analysing the initial flow field; in particular, the velocity at the four VCs is analysed which gives a first hint of where each vortex travels in the beginning. Figure 4 shows the initial positions of the four vortices (marked by asterisks) together with the velocities induced by the surrounding vortices at these positions (indicated by the arrows) for eight different cases. The arrows indicate the direction and strength of this velocity induction; the individual contributions of the vortices V 1 (red), V 4 (green) or the vortex dipole V 2 /V 3 (blue) and their combined effect (black) are depicted. The stream function of an ordinary vortex has circular contour lines and the arrows are tangentials. Hence, the green arrow originating from the red VC, e.g. is perpendicular to the line connecting the red VC with the green one. The contributions of the two IVs (blue arrows), which are counter-rotating, largely cancel out each other far from their centres. Hence, only their combined (and small) effect on the OVs V 1 and V 4 is depicted. On the other hand, the blue arrows that originate from the IVs depict the velocity induction only from the adjacent IV (there is no self-induction). It follows that those blue arrows are perpendicular to the line connecting the two IVs.
The classical single aircraft case is introduced as a reference in panel (a). At both VCs, the velocity is given by u = (0, −w 0 ). With t=0 = 520 m 2 /s and b 0 = 47.3m (as in the formation flight scenarios), w 0 = 1.75m/s follows.
Moreover, the inserted legend introduces the angle φ which defines the propagation direction.
Panel (b) shows the default formation flight simulation. The velocity at the two OVs has a dominant downward component with a small lateral component to the right (black arrows). The contribution of a specific vortex is plotted in the same colour as the asterisk labelling its centre. Note that the blue arrows originating from the OVs V 1 and V 4 show the combined effect of both IVs V 2 and V 3 , whereas the blue arrows originating from V 2 and V 3 depict the velocity induction of the adjacent IV. At each vortex, the combined effect of all surrounding vortices is given by the black arrow. The lengths of the grey arrows inserted in each panel represent a wind speed of 2m/s. We use different scales for the V 1 and V 4 -arrows on the one hand and the V 2 and V 3 -arrows on the other hand reflecting the fact that the IVs move faster than the OVs.
The downward transport is mainly induced by the other OV (green or red arrow, resp.). As noted before, the combined effect of the IVs on the OVs is rather small (blue arrows) and leads to a small counterclockwise rotation of the resulting velocities at V 1 and V 4 . Compared to the single aircraft case, the downward velocities are around a factor of two smaller, as the separation distance between the OVs (given by b 0 + DX ) is roughly twice as large as in the single aircraft case (vortex separation b 0 ). The IVs are transported downward by both OVs.
In sum, the downdraught is twice as large as in the single aircraft case. Yet, the strongest induction comes from the adjacent IV. The IVs mutually induce a strong transport to the right. In total, the IVs travel with more than 5m/s in φ ≈ 110 • direction. The second row shows scenarios with a smaller (panel c)) and larger (panel d)) lateral separation DX compared to the default scenario. We find similar velocity inductions at the OVs. Most notably, the mutual induction of the two IVs changes. In all cases, a strong lateral transport to the right remains. The vertical component, however, might be positive (for DX > b 0 ) or negative (for DX < b 0 ). In total, we find an initial transport of the IVs at a speed of 6.5m/s in φ = 135 • -direction for DX = 45m. For DX = 60m, the vertical velocities cancel out each other and the IVs travel purely horizontally with u = 3m/s. A change in the vertical offset DZ (third row) has similar effects as a DX-variation. Again, the mutual induction of the IVs is affected the most and changes their initial speed and direction. The two remaining cases in the fourth row will be discussed later. We can conclude that in all treated cases the early vortex evolution is characterised by a fast movement of the IVs towards the OV V 4 . In general, the movement is always towards the OV of the lower vortex pair, which is produced by the LAC. Figure 5 shows the axial vorticity field averaged along flight direction at different vortex ages. The positions of the VCs are tracked by applying a simple vortex identification method based on finding local vorticity extrema.
Four simulations with different DX − DZ combinations are selected to demonstrate the diversity and complexity of possible vortex evolutions. The first row shows vorticity fields of the default simulation at t = 16s, 1, 3 and 7.5min. In the beginning, both IVs approach V 4 . Soon, V 3 gets closest to V 4 and those two vortices feature a strong mutual induction. The V 3 &V 4 vortex pair moves quickly away; V 3 propagates more than 300m in the first minute. The V 3 ,V 4 -trajectories bend upwards in counterclockwise fashion. After 1min both vortices are around 50m above flight level, and more than 60m above their creation level. Soon after, they slow down and seem to dissolve, at least our simple method based on longitudinally averaged vorticity evaluation is not able to identify them any longer. However, the 3D representation of vorticity in Fig. 3 has already revealed that the two vortices indeed have dissolved at that stage.
The vortices V 1 V 2 also form a pair and mutually induce a slow downward transport. After 3min, the vertical displacement of V 1 and V 2 is only 70m and 110m, respectively. Both vortices are long-living and, in particular, V 1 can be tracked over 9min (the latest time shown here is 7.5min, though). Whereas V 2 more or less ceases to move and the VC resides around z = −100m, the direction of the V 1 -motion reverses at some point. Then V 1 starts to move upwards and eventually reaches its initial level z = 0m.
The second row shows the simulation where DZ is raised from 14m to 25m. The vorticity distribution after 16s looks similar to the latter case. The relative positions of the vortices V 2 , V 3 and V 4 are only slightly different. Again, V 3 and V 4 build a vortex pair. However, the small initial differences qualitatively change its trajectory and the vortex pair travels sideways. After 1min, V 3 is displaced by nearly 200m in lateral direction. Soon after, the motion slows down. After 4.5 min the vortex pair propagated only another 100m in lateral direction and moved slightly upwards. Eventually the vortices dissolve right above their creation level. As in the latter case, the vortices V 1 &V 2 build a vortex pair and travel straight downward during the first 2min. Contrary, two the latter case, both vortices slow down, remain at around z = −100m and move slightly to the left.
The third row shows the simulation where the original DZ is retained and DX is increased by 10 to 60m. The trajectories are qualitatively similar to case 2. Compared to this case, the final position of the V 3 &V 4 vortex pair is, however, only 200m to the right (instead of 300m in case 2) and around 80m lower. Moreover, the V 1 &V 2 vortex pair travels further down and ends up 150−200m below the creation level.
The forth row shows a case with smaller DX (45m instead of 50m in case 1). Despite this small shift in lateral separation, the vortex evolution differs largely and features new aspects not seen in cases 1 to 3. At t = 16s, the vortex V 2 is in between of V 3 and V 4 . V 2 and V 4 have the same sense of rotation, the weaker vortex (V 2 ) merges in a turbulent way into the strong vortex (V 4 ), and V 2 cannot be tracked beyond t = 30s. Then, V 3 spins around the strong vortex V 4 in counterclockwise sense of rotation. After 70s, the V 3 trajectory completed a full cycle around V 4 and approaches V 1 . V 1 and V 3 have the same sense of rotation. And again, the weaker vortex (V 3 ) merges into the stronger vortex (V 1 ) with the effect that V 3 can be tracked for less than 2min. After 2min, only V 1 and V 4 are present at an altitude of around z = −100m. They continue to descend, though the descent slows down and comes basically to a halt after 3min. The final vertical displacement is less than 150m.
The four simulations exhibit qualitatively different vortex evolutions despite their rather small (< 10m) initial differences. To sum up: What are new aspects compared to the classical vortex descent behind a single aircraft?
We find a strong initial interaction of the two IVs and their lateral propagation (towards the OV of the lower vortex pair generated by the LAC). Afterwards, different phenomena can be observed in the various simulations: a strong lateral transport of the LAC vortex pair, its dissolution far above the flight level, a merging of co-rotating inner and outer vortices, or a reduced descent of the FAC vortex pair. In one case, a single vortex could be tracked over 9min, much longer than in the single aircraft case. Notably, the wake vortex evolution depends very sensitively on the initial lateral and vertical offset.

Plume dimensions
In this section, the spatial distribution of ice crystal number concentrations is analysed. Figure 6 juxtaposes the four FF simulations discussed before and the classical SA (single aircraft) case taken from U2014. For the moment, the ice crystals can be regarded as a passive tracer, i.e. changes in the contrail structure are only due to transport and not due to microphysical loss processes. Hence, we use the terms plume and exhaust instead of contrail and ice crystals in the following. We shortly recall the spatial plume initialisation. It consists of two circular discs with a uniform concentration field where the centres are collocated with the initial LAC VCs. The plumes of the FAC are given by Gaussian distributions. Their centres and peaks, respectively, are aligned with the lateral engine position and are inboard of the FAC VCs. The exhaust gets entrained into the nearby vortices within the first 10s see e.g. Fig. 33.2 in (65,66) . Then, most of the exhaust is trapped inside the vortices and hence the plume expansion is roughly along the vortex trajectories. In the single aircraft case, the vortex descent leads to a substantial vertical plume expansion. In a thermally stable atmosphere, baroclinic instabilities occur around the vortices and exhaust is continuously detrained form the vortex system (67) . The detrained exhaust moves up again due to buoyancy and makes up a curtain between the actual and the initial vortex position (68,69) . After 5min the highest plume concentrations can be found around the original emission altitude. Moreover, the lower part of the plume broadens between t = 3 and 5 min since the vortex tubes start to meander along flight direction (over which is averaged here) as a result of the Crow instability (64) . In this example, the plume is eventually 480m deep and 330m broad. After 5min the vortices have dissolved and the simulation ends (hence, no panel for t = 9 min exists for the single AC case). The further plume expansion is dominated by atmospheric processes and its analysis is not part of this study.
In the FF scenarios, the vortices live longer and the third row shows the panels for t = 9 min. Moreover, all plot axes (including the single AC cases) use the same scale. From the different panel sizes used in the SA case on the one hand and FF cases on the other hand, it becomes directly apparent, that the plume dimensions are qualitatively different. For t = 3 and 5 min, the displayed FF plume structures in Fig. 6 are consistent with the vortex patterns found in Fig. 5. The vertical/lateral expansion is smaller/larger than in the SA case. After the vortices have dissolved and parts of the exhaust moved upwards due to buoyancy (t = 9 min), the plumes are at most 280m deep (cf. to 480m in the SA case) and 500−700m broad (cf. to 330m in the SA case).

Sensitivity to lateral and vertical offset and relative humidity
Sections 3.2 and 3.3 discussed in detail several selected simulations. Contrarily, this section will give a systematic overview over all simulations and summarises physical sensitivities to the lateral offset DX , vertical offset DZ and relative humidity RH i . Whereas DX and DZ strongly affect the vortex evolution (as seen before), RH i affects only ice microphysical processes in the contrail.
The contrail structure is analysed by presenting profiles of ice crystal number in transverse and vertical direction. From this the contrail width and depth can be derived. The contrail depth, in particular, increases very quickly due to the wake vortex descent, much faster than by sedimentation or turbulent diffusion in later stages. Figure 7 shows vertical profiles of ice crystal number for various values of DX (1st column), DZ (2nd column) and RH i (3rd column) for three different instances of times (3, 5 and 7-8min). The discussion of the 4th column will be deferred to the appendix. The vertical   Table 2. The colours are defined in Fig. 9; the black solid curve in each column shows the same default simulation. The RH i -column additionally shows vertical profiles of the classical SA case for the same range of RH i -values (simulation data taken from U2014). The most obvious aspect in the Fig. 7 is the difference between the FF scenarios on the one hand and the SA scenarios on the other hand. As already noted for the exemplary case in Fig. 6, the vortex pair in the single AC cases descends further down than the vortex system in the FF cases. After 3min, the SA plume extends from flight level down to 340m below flight level. In any FF scenario, the plume does not extend beyond 170m below flight level and much more ice crystals are present above flight altitude. After 5min, the SA plumes reach z = −400m and extend only moderately into regions above flight altitude, yielding a contrail vertical extent of 480m in the maximum case. In any FF scenario, no substantial further downward movement of the plumes is apparent and all plumes lie in the region between z = −190 and 170m. The thickest contrail is about 300m thick and several thinner contrails are barely thicker than 200m.
The RH i -column shows a strong RH i -effect on contrail ice crystal number which has been long known for the classical case (31,33) . Downward moving air expands and heats adiabatically. Thereby the saturation vapour pressure increases, the local relative humidity decreases and leads to sublimation and loss of ice crystals. The weaker the ambient supersaturation (i.e. the closer RH i is at to 100%), the more ice crystals get lost during vortex descent. For the classical SA case, all ice crystals in the lower part of the plume vanish for RH i ≤ 110% (red, green, black curves) and the contrail vertical dimension is reduced compared to the moister cases  Table 2. The colours are defined in Fig. 9; the black solid curve in each column shows the same default simulation. RH i = 120% or 140% (blue and brown curve). In the FF scenarios, again fewer ice crystals survive, the lower the ambient RH i is. Yet, the RH i -effect is not as pronounced as in the SA cases and the contrail depth is about the same for all RH i -values. This is due to the fact that the vertical plume displacement is less than in the SA case and the plume adiabatic heating is not as strong. Next, we discuss the sensitivity to DX and DZ (1st and 2nd column). Unlike to RH i , however, no systematic trends can be found in the sense that e.g. contrails become deeper for larger DX or DZ. For example, after 3min the contrail top varies from z = 0 to z = 150m among the four DX-simulations and from z = 50 to z = 150m among the four DZ-simulations. The contrail with the highest top (at z = 150m) occurs for the intermediate values DX = 50m (out of the interval [45, 60m]) and DZ = 7m (out of the interval [0, 25m]). Another example of non-systematic results can be seen in the DZ-panel for t = 5 min. The vertical distributions for DZ = 0 and 25m are nearly the same. Moreover, the DZ = 14m-contrail has a similar vertical extent, yet with somewhat lower ice crystal numbers. The contrail with DZ = 7m differs qualitatively from the three other DZ-cases, as the contrail top is located around 100m higher. The vertical profiles after 7 − 8 min show that in all FF cases, the contrails are around 170 to 250m deep; the variation within the DX-series is even less than 30m. This small spread is remarkable as we have discussed and pinpointed the quite different vortex evolutions for small changes in DX or DZ.
Next, transverse profiles for the same set of simulations are presented in Fig. 8. Again no trends with DX and DZ are apparent and we only report and shortly discuss results for t = 5min and t = 7 − 8min. In the SA setup the single vortex pair has decayed after < 5min and simulations stop around t = 5 min. Hence, the upper row of the figure enables the comparison between the FF and SA scenario, whereas the lower row reports the contrail dimensions after vortex decay for the FF scenarios. At initialisation, the FF contrails are around 130−140m  Table 2.
(depending on DX) and SA contrails are around 90m broad. Within the first minute, a strong lateral transport occurs in some FF setups and the FF contrail width ranges from 200 to 320m (at t = 1 min, not shown). Due to the different strength of the lateral transport in the individual FF cases, the transverse distribution differs strongly from case to case. We find distributions with one or two distinct peaks. The strongest peaks can occur at around x = 0m in one case or at x = 300m in some other case. At t = 5 min, the contrail width varies from 350m for narrow FF contrails to 570m for the broadest FF contrail. A variation of RH i changes the FF contrail width only weakly, whereas for the SA contrails the width ranges from 130 to 290m due to a Crow instability triggered vortex meandering in the lower plume part. Hence, typical SA contrails are usually not as broad as such behind a formation. FF contrails continue to spread after t = 5 min and for t = 7 − 8 min contrail width ranges from 400 to 680m. Integrating the vertical/transverse profiles over the vertical/transverse dimension gives the total ice crystal number N per meter of flight path. Figure 9 shows the temporal evolution of the normalised total ice crystal number f N (t) = N (t)/N (t = 0). Note that in the FF scenarios N (t = 0) is twice as large as in the SA scenarios. The quantity f N diminishes over time due adiabatic heating induced ice crystal loss. The curves flatten out after several minutes indicating that the vortices stopped descending and pushing the local plume RH i below the saturation level. In the chosen classical SA case, this happens after 3 to 4min (the time of vortex break-up depends on the AC type and atmospheric stability, though, see UG2014). In some FF scenarios, the ice crystal loss gets negligible after around 6min. But in many cases crystal loss does not fully vanish and occurs over the complete simulation period, though at a reduced rate towards the end. Again no trends with DX and DZ can be observed and the fraction of surviving ice crystals f N ,s lies between 0.4 and 0.6 (surviving here refers to the f N -value at the end of the simulation). RH i is the dominant parameter for f N ,s . Generally, a larger fraction of ice crystals survives in the FF than in the SA scenarios for a given RH i -value. When RH i is reduced from 140% to 120% and further down to 110%, f N ,s strongly decreases in the SA case (from 0.92 to 0.67 and further down to 0.28). In the FF scenarios, on the other hand, a reduction of RH i = 140% down to 120% (the supersaturation RH i − 100% is actually halved) has only a moderate effect and f N ,s goes down from 0.96 to 0.89. For both RH i -values, the vortex descent does not suffice to produce a substantial crystal loss. Only, a further reduction down to 110% makes the atmospheric buffer effect so small that the reduced vortex descent of FF cases produces local subsaturations that are sufficient to lead to a substantial loss of ice crystals (f N ,s = 0.62).

Robustness of findings
Various further experiments and a model comparison have been carried out in order to check the robustness of the above findings. These efforts are shortly summarised here and a full description is deferred to the appendix. In a grid sensitivity study, where the resolution and domain size along flight direction was varied, we find that the grid choices have a minor impact on the quantities of interest. Moreover, we perform four different realisations of the same simulation by simply shifting the turbulent background fields in lateral direction. We find small turbulence induced uncertainties implying the significance of the physically induced differences seen in the results section. Section 2.2.2 described various versions of IV initialisations. The initialisation choices have a non-negligible impact on the simulation results. Reducing this uncertainty will be an important future task and the way forward will be laid out in the discussion section. Finally, MGLET simulations are repeated for different DX-DZ combinations. In all three MGLET simulations we find patterns consistent with those of the EULAG simulations. Keeping in mind the diversity and complexity of the observed vortex evolution for small variations of DX or DZ, it is very convincing that the patterns are similar in both models and this boosts the confidence in the robustness of the simulation results.

DISCUSSION
In this section we discuss implications of the presented results.
Formation flight would probably reduce the frequency of potentially hazardous wake encounters at cruise levels for several reasons. The frequency of wake-vortex encounters increases with the square of air-traffic density (70) . Organising the sky by routinely matching up several aircraft in a formation the flight routes of individual aircraft are not randomly distributed any longer and the effective air traffic density is reduced. Second, the present study demonstrated that, at least for two aircraft configurations, the descent speed and the final vertical displacement of the vortices in a FF scenario are smaller than in the classical SA scenario. Hence, the probability of interferences across different flight levels (the main reason for in-cruise wake vortex incidents (70)(71)(72) ) is likely to be reduced. Moreover, the experienced roll moments behind an AC formation wake may differ from those behind a single aircraft wake and may impact the severity of encounters. The latter aspect could be elaborated in a follow-up study based on the simulated flow fields. One aspect that may counteract those three effects is the fact that wake vortices appear to be more long living than in the single aircraft case.
Next, we discuss the implications of FF on long-living contrail-cirrus, which may persist for many hours (26) and have a substantial climate impact on global scale (23) . Regarding the long-term evolution over many hours, where sedimentation becomes important, the vertical expansion during the first several minutes can be viewed as a 'sudden' event. The spreading rate, i.e. the rate with which the contrail width increases with time, is proportional to vertical wind shear (which is an ubiquitous phenomenon in the upper troposphere) and to contrail depth. Moreover, sedimentation can be a limiting factor of contrail life time (73) . The fewer ice crystals are present initially in the contrail, the bigger they grow during contrail aging and the faster they fall out. Moreover, fewer, but larger ice crystals (with the same total ice mass) lead to an reduced extinction of shortwave radiation. For those reasons, the contrail vertical extent and the total ice crystal number after vortex decay are important quantities that can affect the contrail-cirrus properties over hours (25,30) . UG2014, e.g. simulates contrails produced by various aircraft types and initial differences in ice crystal number and contrail depth lead to quantitative differences between the contrail-cirrus properties that remain over the total simulation period of 6h.
Earlier studies (29) showed substantial saturation effects, when contrails were produced in close proximity and form a contrail-cluster. Compared to isolated contrails, the overall effect on radiation by an individual contrail in a cluster is smaller. Those as such contrails compete with neighbouring contrails for the available water vapour and hence their deposition growth is limited. In the formations studied here, two aircraft produce one 'multi'-contrail. This case can be interpreted as the limiting case of the aforementioned contrail cluster simulation with initial contrail separation distances between 5 and 20km, as the distance of the two involved contrails tends to be much smaller here. So small even, that the distance is not well-defined any longer, as the contrails start to interact from the beginning of their existence. Hence, one can expect strong saturation effects in the case of formation flight simply due to the fact, that a single multi-contrail is produced by two AC. Without formation flight two independently produced contrails most likely occupy different air masses and can evolve undisturbed from each other.
The extent of saturation may be also affected by the fact, that a FF contrail contains more than twice as many ice crystals and is shallower than a single aircraft contrail. Contrail-cirrus simulations based on young FF and SA contrails will be presented in a follow-up study that investigates in more detail the climate benefit of formation flight due contrail modification.
As a last point of the discussion, a road map for reducing uncertainties associated with the inner vortex initialisation is outlined. In this study, temporal LES have been performed both with EULAG-LCM and MGLET, where a homogeneous vortex flow field together with periodic boundary conditions along flight direction is used. MGLET also features more advanced initialisation features, in particular for near ground wake vortex applications. A one-way coupling of RANS and LES simulations provides a methodology of simulating an aircraft flight through a computational domain, generating a wake in a realistic way (74,75) . For this purpose a pre-computed high-fidelity steady RANS flow field is swept through the LES domain. This is done to initialise the aircraft wake in the LES domain and is referred to as 'spatial LES'. This method includes all stages of wake-vortex evolution, from roll-up to vortex decay and can capture, e.g. the touchdown manoeuvre and the wing in ground effect at varying vortex generation heights above the ground. The results have revealed valuable insights for cruise flight (75) as well as flight in ground proximity during the landing phase (76) . However, the problems addressed in this work cannot be tackled with a one-way coupling, since it requires the following aircraft response on the wake of the leading aircraft. Therefore it is not applied here. Recently, a two-way RANS-LES coupling was developed, where a RANS code, simulating the aircraft flow, is coupled bidirectionally with the LES code MGLET. This allows simulating a realistic aircraft flight through a turbulent environment in a ground fixed domain. Having a leading and a following aircraft in that framework will provide a more realistic approach to formation flight, in particular valuable insight on the near field wake-vortex roll-up and early interaction can be obtained.

CONCLUSIONS
• For the first time, large-eddy simulations of the far field wake vortex evolution behind an aircraft formation have been performed, hereby focusing on a stratified atmosphere typical of cruise altitudes.
• The two LES models EULAG-LCM (43,53) and MGLET (77) have been employed separately for such simulations. Both finite-difference codes have been used for wake vortex modelling in the recent past e.g. (1,2,36,38) . Moreover, the LES model core EULAG is fully coupled with the Lagrangian ice microphysics code LCM, which allows simulating the contrail evolution; in a reduced version LCM can be used for the Lagrangian transport of a passive tracer (35) . • Besides evaluating vorticity-related quantities and tracking of vortex centres, the present study investigates the aircraft plume dispersion which depends on and also reveals the wake vortex movement. their descent leading to a substantial plume/contrail expansion in vertical direction. In formation flight scenarios, wake vortex behaviour was found to be much more complex and also diverse. In the beginning, the two inner vortices always propagate towards the outer vortex of the lead aircraft. The exact propagation direction depends on the initial relative position of the inner vortices (i.e. on DX and DZ) and determines the further fate. Once the two inner vortices get close to the outer vortex, remarkably different things occur in the various simulations. In one simulation setup, it happened that the co-rotating inner vortex soon merges into the outer vortex. Later on, the two other vortices, both with the opposite sense of rotation, also merge. In several other simulations, on the other hand, the outer vortex forms a vortex pair with the counter-rotating inner vortex. Subsequently, this vortex pair tends to travel sideways or upwards. The two other vortices most often form a slowly descending long-living vortex pair which could be tracked for more than 5min. In one case, a vortex tube is still apparent after 7.5min and a single vortex can even be tracked over 9min. The wake vortices behind a single aircraft would have been fully dissolved after 5min in the same atmospheric environment. • Notably, CO 2 emissions are only one contribution to the aviation climate impact among several others like contrails and emission of water vapour and nitrogen oxides (20,21) , where contrails probably have the largest radiative forcing (22,23) . All those components could be affected by formation flight. The (classical) early contrail evolution is characterised by vertical expansion due to vortex descent and by ice crystal loss. The latter is caused by adiabatic heating in the downward sinking primary wake, which leads to a local reduction of relative humidity and ice crystal sublimation. • With EULAG-LCM the early contrail evolution and its sensitivity to relative humidity over ice RH i was investigated. We find that, compared to the classical case, fewer ice crystals are lost as the adiabatic heating is not as pronounced. Hence, the formation flight contrail contains more than twice as many ice crystals as the single aircraft contrail. Moreover, the contrails are shallower and broader which are a direct consequence of the wake vortex movements. • A follow-up numerical study that uses the present simulation data as input will focus on the contrail-cirrus evolution over several hours and assess the climate benefits by contrail saturation effects in formation flight scenarios.

NUMERICAL CONSISTENCY
Whereas the results section focused on the impact of physical parameters, the appendix presents mostly numerical sensitivity tests and a model comparison in order to demonstrate the robustness of the presented results. The first subsection discusses the uncertainties associated with the uncertain initialisation of the two IVs. The second subsection investigates the sensitivity to grid parameters for the EULAG default simulation, whereas the third subsection highlights irreducible uncertainties associated with turbulence. The appendix concludes with a comparison of EULAG and MGLET results for three selected cases.

Uncertainties from the inner vortex initialisation
In Section 2.2.2 we described various versions of IV initialisations. The simulation results are presented in the right-most columns of Figs 7, 8 and 9. In the default case, the two IVs are damped in R weak = 15m-circles around the VC and turbulence is enhanced in a box around those two vortices (see yellow box in Fig.1). This simulation is shown in black (labelled by '15' in the legend).
In a first sensitivity test, we increased R weak to 25m (label '25'). The bottom right panel of Fig. 4 illustrates how an increase of R weak affects the initial flow field and propagation direction of each vortex. The induced velocities at the OVs are unaffected as their distance from the IVs is larger than 25m and IVs are not damped at that radii, irrespective of R weak being 15 or 25m. The mutual induction of the IV is, however, reduced (see blue arrows). Together with the unchanged contribution of the OVs on the IVs, the propagation direction of the IV has a stronger downward and a weaker horizontal component (see black arrows).
We further tested two variants, where the vortex strength was reduced for all radii (i.e | 2 | and | 3 | are reduced), and not only locally close to the core as before. In one test, both vortices are similarly weakened ( 2 = − 3 = 0.5 t=0 , label '5050'). In a second test, the reduction is unsymmetrical ( 2 = 0.6 t=0 and 3 = −0.4 t=0 , label '6040'). The bottom left panel of Fig. 4 shows again the immediate consequences on the initial vortex propagation. Unlike to the R weak -change before, the full damping assumption also changes the propagation of the OVs. The IVs pull the OVs to right. If IV circulations are lowered, this effect is reduced and the OV propagation tends to become nearly purely vertical. The mutual induction of the IVs and their propagation is similar to the R weak = 25m-case.
Furthermore in yet another test (label 'ET0'), turbulence was not enhanced. And last, we ran a simulation with the overly simplistic case with no damping at all, i.e. four identical vortices in the beginning (label '00').
A simple observation from Figs 7, 8 and 9 is that the inclusion/omission of enhanced turbulence has a negligible impact on the vortex evolution, the spatial distribution of the ice crystals and ice crystal loss. Hence, the ET0-simulation will not be discussed any further in the following.
We recall the VC evolution of the default case as shown in Figs. 3 and 5. The IVs V 2 and V 3 move towards the right OV V 4 . V 4 picks up V 3 and those two vortices form a vortex pair that decays quickly. Their trajectories are next to each other and bend upwards in a counterclockwise fashion. V 1 and V 2 also form a vortex pair that is slowly moving down and is more persistent than the other pair. In the various IV sensitivity studies, qualitative deviations from this rough pattern occur (not shown). In one simulation, e.g. the vortex V 2 remains longer under influence of V 4 and does not form a pair with V 1 . In another simulation, the co-rotating vortex V 2 merges into V 4 at some point, the vortices V 1 and V 4 form a persistent and slowly decaying vortex pair. The implications of those differences can be seen in the transverse and vertical ice crystal number profiles. The vertical profiles across the IV simulation series show some variation. The transverse profiles show a small spread (neglecting the out-lier simulation '00') and non-zero ice crystal numbers are confined to the area x = −200m to 250m. Also the extent of ice crystal loss depends on the IV initialisation. In the default case, around 60% of the ice crystals survive, whereas in most other IV simulations around 45% of the ice crystals survive. In general, it seems that the spread in the profiles is smaller across the IV series than across the DX -or DZ-series. Yet, the initialisation choices have a non-negligible impact on the simulation results. Reducing this uncertainty will be an important future task as detailed in the discussion section.

Grid sensitivity
This section presents a grid sensitivity study for the default simulation (DX = 50, DZ = 14m, RH i = 110%). First, the vertical and lateral extent have to be chosen large enough to avoid interactions across the periodic lateral boundaries or with the rigid boundaries on the top/bottom. This was assured by choosing a fairly large domain with L x = 768m and L z = 600m compared to previous classical wake vortex studies (36) and similar to a study accounting for cross-wind effects (78) . As in previous studies the grid boxes use x = z = 1m. Along flight direction, the domain in classical wake vortex simulations is chosen such that one Crow modes can develop, i.e. around 400m for an AC with a wing span of 60m. Due to the strong interaction of the two IVs one cannot expect that a Crow mode with 400m long oscillations will occur in our FF simulations. Hence, we tested domains with various numbers of grid points n y = 396, 198, 99 and 66 and grid sizes y = 2, 1 and 0.5m. The section grid parameters in Table 2 and the legend in Fig. 10 lists all six combinations of L y = n y · y and y that were tested. The largest simulation domain has a length of 792m and may contain a Crow oscillation of the OVs which have distance of around 2b 0 or several Crow modes of two other vortices with a smaller separation.
Concerning the transverse distribution of ice crystal number, shown in the bottom row of Fig. 10, we find an excellent agreement among all six simulations. Also the vertical profiles shown in the top row match well in a qualitative sense, yet with increasing time the distributions of the various simulations differ moderately close to the contrail top. Moreover, no distinct oscillations along flight directions could be identified by inspecting longitudinal profiles (not shown). This finding is corroborated by the 3D vorticity fields in Fig. 3.  Table 2 (section grid parameters)). Figure 11. Temporal evolution of normalised ice crystal number. Left: Grid sensitivity experiment, simulations and colours as in Fig. 10. Note that the cyan curve is hidden by the red curve. Right: Turbulent realisation experiment, simulations and colours as in Fig. 12.
The negligible impact of the grid on the crystal loss mechanism and the total ice crystal evolution is exemplified in Fig. 11. All in all, we find that the choice of the grid has a minor effect on the quantities discussed in the results section. Hence, we decided to use the most inexpensive setting with L y = 132m and y = 2m which helped to keep the computational requirements in an acceptable range and allowed to perform several physical and numerical sensitivity studies.

Turbulent realisation
In this sensitivity test, we use the default grid and also keep the default values of all other parameters unchanged. We perform four different realisations of the same simulation by simply shifting the turbulent background fields in lateral direction. Although the statistical turbulence properties are unchanged, the vortex system may evolve differently due to slightly different local interactions of the vortices and turbulence. The spread of the four simulations demonstrates the irreducible uncertainties associated with turbulence. If the simulation results of two simulations with different parameter settings (as presented in the results section) differ, those differences are only significant if they are larger than the spread given by the various realisations of a specific simulation. We find small turbulence induced uncertainties ( Fig. 12 and right panel of Fig. 11) implying the significance of the physically induced differences seen in the results section.

Comparison of EULAG-LCM and MGLET simulations
This section presents a comparison of EULAG-LCM and MGLET simulations. The LES model MGLET, developed at the Technical University of Munich (77) , solves the incompressible NS equations. The kinematic viscosity is given as the sum of molecular viscosity ν and eddy viscosity ν t , determined by means of a Lagrangian dynamic subgrid scale model (79) . The momentum equation and the continuity equation are solved by a finite-volume approach, using a fourth-order compact scheme by Hokpunna and Manhart (80) . A split-interface algorithm is used for the parallel solution of the tri-diagonal system resulting from the compact scheme (81) . A third-order explicit Runge-Kutta method is used for the time integration. The simulations are performed in parallel, using a domain decomposition approach. MGLET uses the same turbulent background fields at initialisation as EULAG and superimposes the vortex fields in an identical way. EULAG uses Lagrangian particles with discrete positions for the representation of ice crystals and their advection. MGLET is not coupled to a ice microphysics module, but it offers the advection of a passive tracer represented by Eulerian field. Unlike to Lagrangian methods, Eulerian advection schemes do not perform well when strong gradients appear in the tracer distribution. Hence, in MGLET the plume initialisation with circles of uniform concentrations uses some smoothing in a transition layer. Due to operational constraints, MGLET simulations use a domain width of L x = 600m instead of 768m and the simulated time is 5min. Three of the four scenarios discussed in detail in Figs 5 and 6, namely the default simulation, the DZ = 25m-simulation and the DX = 45msimulation, have been repeated with MGLET. Figure 13 depicts the vorticity distribution and tracks the movement of the VCs in the new simulations. The plot layout is analogous to the EULAG counterpart in Fig. 5, except that one simulation and the last point of time are discarded. In all three MGLET simulations we find patterns consistent with those of the EULAG simulations. In the default case, the V 3 ,V 4 -trajectories bend upwards in counterclockwise fashion and the vortices V 1 &V 2 form a pair with a slow downward transport as already seen in EULAG. In the DZ = 25m-case, we again find a strong purely horizontal transport of the V 3 &V 4 vortex pair. And in the DX = 45m-case, we observe a merging of the co-rotating vortices V 3 and V 4 as in EULAG. Figure 14 overlays the vortex positions (left column) of EULAG (solid) and MGLET (dotted) and evaluate the distances between the EULAG and MGLET vortex positions (right column). The differences between the two models increase over time, but typically remain smaller than 50m after 5min for all four vortices. Keeping in mind the diversity and complexity of the observed vortex evolution for small variations of DX or DZ, it is very convincing that the patterns are similar in both models and this boosts the confidence in the robustness of the simulation results. Figure 15 shows the MGLET passive tracer distribution similarly to the EULAG ice crystal number concentrations in Fig. 6. Since the vortex trajectories agreed well, it is not surprising that also the plume cross sections look similar. Again, the plume width and vertical extent are inserted in each panel. The given values mostly deviate less than 50m from EULAG values. Considering the different advection treatment, the agreement is excellent. Finally, Fig. 17 juxtaposes transverse and vertical profiles of MGLET (right) and EULAG (left). As default, RH i = 110% is used in the EULAG simulations. As a consequence, ice crystal loss occurs in all simulations of the DX and DZ sensitivity series and in particular in the three simulations