The Contrail Mitigation Potential of Aircraft Formation Flight Derived from High-Resolution Simulations

Formation flight is one potential measure to increase the efficiency of aviation. Flying in the upwash region of an aircraft’s wake vortex field is aerodynamically advantageous. It saves fuel and concomitantly reduces the carbon foot print. However, CO2 emissions are only one contribution to the aviation climate impact among several others (contrails, emission of H2O and NOx). In this study, we employ an established large eddy simulation model with a fully coupled particle-based ice microphysics code and simulate the evolution of contrails that were produced behind formations of two aircraft. For a large set of atmospheric scenarios, these contrails are compared to contrails behind single aircraft. In general, contrails grow and spread by the uptake of atmospheric water vapour. When contrails are produced in close proximity (as in the formation scenario), they compete for the available water vapour and mutually inhibit their growth. The simulations demonstrate that the contrail ice mass and total extinction behind a two-aircraft formation are substantially smaller than for a corresponding case with two separate aircraft and contrails. Hence, this first study suggests that establishing formation flight may strongly reduce the contrail climate effect.


Introduction
Formation flight (FF) is a well-known strategy of migratory birds in order to improve their aerodynamic efficiency, save energy and increase their range [1][2][3]. Similarly, FF can increase the performance in the civil and military aviation sector. In addition to close FF (with separations of a few wing spans in flight direction) in the military sector, extended FF (with separations of 10 to 40 wing spans) is a viable option for the commercial sector. Follower aircraft (AC) encounter uplift from flying in the upwash region created outboard of a leading AC. Numerous numerical, wind tunnel and real flight studies, e.g., [4][5][6][7][8][9][10][11][12] demonstrate that the induced drag is reduced, the lift-to-drag ratio increases and fuel consumption is lower. One goal of these studies is to find the lateral and vertical offset of the AC positions with maximum benefit (the so-called sweet spot).
Fuel savings of around 10% can be expected during such formations (on average over all participating AC, not only for the follower AC; see exhaustive list of references above). In order to establish formations in the airspace, re-routing is required. Clearly, re-routing induced fuel penalties should be substantially smaller than the fuel savings during the actual formation. Xu et al. [13] find net fuel burn reductions of nearly 8% and 6% for different network sizes of cooperating airlines.
Fuel savings directly translate into lower CO 2 emissions and would reduce the climate impact of aviation. Besides the emission of carbon dioxide, the formation of contrails and the emission of nitrogen oxides cause a substantial aviation radiative forcing [14,15] which would be both affected by the introduction of FF. The contrail radiative forcing (RF) is probably larger than the RF of the total accumulated CO 2 emissions from aviation [16,17] and the present study focuses on the contrail mitigation potential by FF.
Contrails grow and spread by the uptake of atmospheric water vapour and the initial emission of water vapour contributes a minor fraction to the ice mass of the aged contrail-cirrus [18] and also to its radiative effect. Saturation effects can be expected when contrails of two or more aircraft are produced in close proximity as they compete for the available atmospheric water vapour and mutually inhibit their growth. Hence, their overall effect is smaller than in a situation with similar atmospheric conditions where each contrail evolved undisturbed of the others [19]. Note the difference to the CO 2 -effect, where the impact is solely determined by the initial emissions as they simply accumulate.
In a recent study, the early evolution of contrails produced by two AC in formation was analysed and compared to classical contrails behind a single aircraft. Classical single aircraft (SA) contrails grow mainly in vertical direction due to the wake vortex descent in the first few minutes, e.g., [20][21][22]. In FF scenarios the individual contrails merge quickly after their formation into a single contrail (from now on abbreviated as FF contrail). Unterstrasser and Stephan [23] found the wake vortex dynamics to be more complex and diverse as vortex pairs happen to also move upwards and sideways. After vortex break-up, FF contrails are thus not as deep, yet they are broader than SA contrails. Earlier studies already demonstrated that differences in the early contrail properties triggered by differences in wake vortex characteristics can have a long lasting mark. Unterstrasser and Görsch [24], e.g., simulated 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 6 h.
We employ the large-eddy simulation (LES) model EULAG-LCM, which is an established code for performing high-resolution simulations of contrails. In this study, we will juxtapose the evolution of contrails generated by a single aircraft, on the one hand, and by a two-aircraft formation, on the other hand. First, an exemplary simulation is presented, then the contrail-cirrus evolution is discussed for selected atmospheric scenarios. Finally, the study evaluates the extent of saturation in the formation flight scenarios (as introduced above) for a large set of prescribed atmospheric scenarios.

Methods
This section introduces the employed model with its numerical set-up and defines properties that will be used in the later analysis. The setup of the contrail-cirrus simulations, also known as dispersion phase simulations, presented here is in many aspects similar to the one in Unterstrasser et al. [25,26]. Hence, we give only a short description here. More details on the numerical setup are given in the two latter references.

Model
For the numerical simulations, the model EULAG-LCM has been used. EULAG [27,28] is a non-hydrostatic anelastic LES model, which employs the positive definite advection scheme MPDATA [29,30]. The ice microphysical module LCM, based on Lagrangian tracking of ice crystals [31], is fully coupled to EULAG. The model version EULAG-LCM has been used for the simulation of natural cirrus and contrails. Application examples are studies of a mid-latitude cirrus cloud system with a special focus on aggregation [32] and the contrail evolution during the vortex phase [24,33]. Moreover, contrail-cirrus simulations have been presented in Unterstrasser et al. [25] and its interaction with surrounding natural cirrus has been analysed in Unterstrasser et al. [26].
In LCM, ice crystals (ICs) are represented by Lagrangian simulation particles (SIPs), where every SIP represents a large number of real IC with identical properties. In order to reduce the complexity of the present simulations and simplify their interpretation, several LCM components are switched off (like heterogeneous nucleation, aggregation, and radiation). Although the switched-off processes can strongly alter the evolution of natural and contrail-cirrus, e.g., [32,[34][35][36][37], we will give arguments in the discussion section why we do not believe that the omission of those processes introduces systematic biases in the comparison of formation flight and single aircraft scenarios. The basic microphysical processes considered in the present set-up are deposition/sublimation and sedimentation. Moreover, natural cirrus can form by homogeneous nucleation. For the ice crystal habit, we assume hexagonal columns. A (synoptic scale) spatially homogeneous updraught motion is prescribed via an external forcing term in the temperature equation in order to accommodate for the adiabatic temperature reduction. The temperature reduction results in an increase of the background relative humidity and of the excess moisture. Details of the forcing implementation can be found in Unterstrasser and Gierens [36].

Numerical Set-Up
The simulation set-up for the present study is similar to that of previous EULAG-LCM contrail-cirrus simulations [24][25][26]. A two dimensional model, whose domain is perpendicular to the direction of flight and represents some portion of the UT/LS (upper tropospheric/lower stratospheric) region, is used. In the vertical (z-)direction, the domain dimension is 2.5 km. In the horizontal (x-)direction, the domain dimension is 40 km or 80 km (depending on the strength of the vertical wind shear; higher wind shear leads to broader contrails). Uniform grid boxes with sizes dx = dz = 10 m span a regular Cartesian mesh. The total simulated time is 8 h. The dynamical time step ∆t dyn is 2 s or 1.25 s depending on the vertical wind shear and the nucleation time step ∆t nuc is always 0.4 s (nucleation is turned on only in simulations with ∆t dyn is 2 s). Table 1 summarises the default values of fundamental numerical and atmospheric parameters.
The vertical temperature profiles corresponds to a stably stratified atmosphere; we prescribe a Brunt-Väisälä frequency N BV of 10 −2 s −1 , a value typical of the upper troposphere. Background turbulent velocity fields are taken from a-priori simulations and have a root mean square (rms) valuê Even without any turbulence forcing mechanism,û does not drop below 0.1 m s −1 and the turbulence intensity is quasi-constant over the simulated period see Figure 3 in [25]. An ice-supersaturated (ISS) layer of d ISSR ≈ 1.2 km thickness and an initial relative humidity (with respect to ice) of RH i,0 * is prescribed between z = 1000 m and 2000 m (see Figure 1 top).
Above and below this layer, RH i,0 drops to 20% inside 500 m thick transition zones. Below this transition zone, a 500 m thick layer with RH i,0 = 20% follows. The initial relative humidity RH i,0 * is either 110% or 120%. The superscript "*" refers to the constant RH i value in the ISS layer. (Bottom): Temporal evolution of the background relative humidity RH i * (t) at z = 2000 m for various updraught speeds w syn = 1 cm s −1 (brown), 2 cm s −1 (blue) and 5 cm s −1 (red and green). The initial RH i,0 * is either 110% or 120%. For RH i,0 * = 110% only one scenario with w syn = 5 cm s −1 is displayed (red). The solid and dotted lines show scenarios with an adiabatic cooling of 2 K and 4 K, respectively.
The flight altitude of the contrail producing aircraft is at z CA = 1500 m in the middle of the ISS layer resulting in d up = d down ≈ 600 m thick fractions of the ISS layer above and below cruise altitude.
Three different synoptic scale updraught scenarios with w syn either = 1 cm s −1 , 2 cm s −1 or 5 cm s −1 are prescribed. By adjusting the updraught period, the prescribed final adiabatic cooling is either 2 K or 4 K and corresponds to an uplift of roughly 200 m or 400 m, respectively. Table 2 summarises the updraught velocities w syn and durations t updr of the various scenarios. The temporal evolution of the background relative humidity RH i * (t) at z = 2000 m is shown in Figure 1 bottom.
Three cases with RH i,0 * = 120% and one case with RH i,0 * = 110% are selected for display. In the 2 K scenarios, the peak RH i is around 150% or 140% depending on RH i,0 * . This is below or close to the threshold humidity where homogeneous nucleation starts. We deliberately switch off homogeneous nucleation in the 2 K scenarios (otherwise ICs might form in parts of the domain where turbulent fluctuations create RH i values above the threshold) and hence contrails are the only cloud type in those simulations. In the 4 K scenarios, the ambient relative humidity rises well beyond the nucleation threshold and natural cirrus is allowed to form around the contrail. Homogeneous nucleation is initiated preferentially at the top of the supersaturated layer where the nucleation threshold humidity of RH crit ≈ 155% is surpassed first. The onset of cirrus formation t nuc depends on w syn (and to a lesser degree on RH i,0 * ) and is listed in Table 2. During the cirrus formation stage, up to 10 × 10 6 SIPs are generated to resolve the highly non-linear ice nucleation process. The stochastic nucleation implementation and an SIP merging technique as described in Unterstrasser and Sölch [38] are employed to increase numerical efficiency. In the very slow updraught scenario with w syn = 1 cm s −1 natural cirrus would form late in the simulated period and would not affect much the contrail evolution. Hence, 4 K cooling cases are only run for w syn = 2 cm s −1 or 5 cm s −1 .
Finally, we summarize aspects that differ from the original setup in Unterstrasser et al. [25]. In the original setup, a dry layer was also included above the upper transition zone. We found this layer to be irrelevant for the simulation results. Leaving it out, we could reduce the domain height from 3 km to 2.5 km. Originally, the flight altitude was close to the top of the ISS layer at z CA = 2000 m with d up = 200 m and d down = 1000 m. The present setup with d up = d down ≈ 600 m was tested in a sensitivity series called ISS_up. Moreover, it was necessary to slightly adapt the number of horizontal grid points nx from 4096 to 4032 due the changes in the supercomputing environment.

Contrail Initialisation
As in previous studies, the contrail initialisation is based on 3D data of contrail vortex phase simulations. The initialisation of a classical SA contrail is based on input data of Unterstrasser [33], whereas the corresponding simulations of an FF contrail-cirrus produced by a two-aircraft formation start with input data of Unterstrasser and Stephan [23].
The following presentation of the initial contrail properties is rather detailed as we will later see that the extent of saturation turns out to be sensitive to the contrail initialisation.
Generally, the contrail vortex phase features a vertical expansion of the contrail due to the downward moving wake vortex and potentially substantial ice crystal loss due to adiabatic heating. A summary of these processes and their sensitivity to atmospheric and aircraft parameters is given in Unterstrasser [22] and Paoli and Shariff [39].
Unterstrasser and Stephan [23] compare young contrails after vortex breakup produced behind a formation to those of the classical case. Qualitative differences were found: Behind a formation, the wake vortices of both aircraft interact. Often this leads to a strong lateral transport of one vortex pair and moderate sinking of the second pair. Hence, contrails behind a formation are broader, yet their vertical extent is smaller than in the classical case. Moreover, the ice crystal loss is not as pronounced as in the classical case.
All vortex phase simulations used as input in the present study use an ice crystal 'emission' index of 2.8 × 10 14 (kg fuel) −1 and the aircraft characteristics of an A350/B777 aircraft. Shortly after contrail formation, the SA contrail then consists of N 00 = 3.4 × 10 12 ICs per meter of flight path, and in an FF contrail the number is just double as high (N 00 = 6.8 × 10 12 m −1 ). Due to crystal loss, the actual number decreases during the vortex phase. Columns 4 and 3 of Table 3 list the survival fraction f Ns and the resulting IC number N 0 with which the present simulations start. Clearly, N 0 is smaller for smaller RH i,0 * . Moreover, the ice crystal number behind a formation is more than double as high as in the classical case, as fewer ICs are lost during the vortex phase. The factors N 0,FF /N 0,SA are 4.5 and 2.6 for RH i,0 * = 110% and 120%, respectively. Table 3. Initial contrail properties: f Ns is the fraction of ice crystals that survive the vortex phase and is determined by f Ns = N 0 /N 00 . In the single aircraft case N 00 = 3.4 × 10 12 m −1 and in the formation flight scenario N 00 = 6.8 × 10 12 m −1 .  Figure 2 shows contrail IC number profiles in vertical (left) and transverse (right) direction. From the contrail profiles the contrail width W and vertical extent H are inferred and the values are also listed in Table 3 (columns 1 and 2). The solid curves shows the SA cases, where the contrail can extend up to 400 m below the flight altitude (in this plot z = 0 is identified with the flight altitude). For RH i,0 * = 110%, the contrail is considerably thinner since basically all ICs that were trapped in the downward sinking vortex system sublimate completely. In the FF scenarios, the contrail vertical extent ranges between 250 m and 300 m and is substantially smaller than for fully grown classical contrails.

Contrail Parameters
On the other hand, FF contrails are broader than SA contrails (400-600 m vs. 150-250 m). Moreover, large parts of the FF contrail can lie above cruise altitude. For this reason, the contrail is chosen to be located in the middle of the ISS layer in the contrail-cirrus simulations. If we had used the standard setup of Unterstrasser et al. [25] with d up = 200 m, the FF would penetrate into the transition zone, which is not desirable for a meaningful comparison between SA-and FF cases. Moreover, the various red lines illustrate the FF contrail sensitivity to parameters of the formation geometry, namely the lateral offset of the two AC in formation. The number concentrations in young contrails are usually so high that the ambient supersaturation quickly relaxes to saturation as confirmed by contrail observations [40]. Hence, the ice water content scales linearly with the contrail volume and the ambient supersaturation. The listed initial ice mass values I 0 reflect this fact. In particular, the FF cases have a higher initial ice mass as their volume is larger than that of classical contrails.
Next, we will shortly describe the technical procedure of how to incorporate the vortex phase simulation data in the model domain. The Eulerian 3D data (e.g., velocity (u, w), perturbations of water vapour concentration q v and potential temperature θ) are averaged along flight direction and interpolated on the coarser grid and embedded into the enlarged 2D model domain. Lagrangian SIP data contain the information of the ICs and SIPs with similar positions (neglecting the coordinate y along flight direction) and IC sizes are merged. This reduces their overall number from several tens of millions in the 3D simulations to around 1.2 × 10 6 in the present approach. Unterstrasser and Sölch [38] showed that contrail-cirrus is sufficiently well represented by 1.2 × 10 6 SIPs for the type of analysis carried out in this study. Figure 3 exemplarily displays an initial 2D field, where the black box encompasses the area in which the 3D vortex phase data was inserted. In this example, relative humidity is shown. The RH i values in the contrail are apparently lower than in the supersaturated environment, where the local RH i values fluctuate around RH i,0 * = 120%.  The SA simulations start with a 5 min old contrail as in previous studies. Behind a formation the wake vortices live longer and the present simulations start with a 7 min old FF contrail. Sensitivity tests using a 5 min old FF contrail showed similar results implying that the FF contrails basically reached their characteristic shape after 5 min and the IC loss has come to an end.
In the literature, the terms contrail and contrail-cirrus attempt to distinguish linear contrails from aged contrails that lost their line shape and whose appearance resembles those of natural cirrus. However, no clear definition can be made on what should be called contrail and what contrail-cirrus. In the following, the usage of the terms contrail and contrail-cirrus is not physically motivated; it should simply ease the description of vortex phase simulations to simulate contrails, whereas dispersion phase simulations simulate contrail-cirrus (CC).

Set of Simulations
This subsection gives an overview over the complete set of simulations and clarifies which parameters and combinations of them are varied. For any given atmospheric scenario, simulations with SA contrail initialisations (REF) and FF contrail initialisation (FF) are performed.
As noted in the previous subsection, there are qualitative changes in contrail dimension and IC number between the REF and the FF initialisation. The later CC properties are affected by those initial differences as will be demonstrated in Section 3. To disentangle the effects of a variation of contrail dimension, on the one hand and of a variation of IC number, on the other hand, a third type of simulation referred to as NNN is introduced. The NNN simulations are based on the SA contrail initialisation and the contrail structure is thus identical to the REF case. The IC number concentrations are uniformly upscaled, such that the total IC number matches that of the corresponding FF scenarios. This upscale factor is 4.5 or 2.6 for RH i,0 * = 110% or 120%, respectively.
The CC evolution and properties depend on many atmospheric parameters; several of them are chosen to determine what we call the atmospheric scenario. These are the initial relative humidity RH i,0 * , the final adiabatic cooling ∆T, the updraught speed w syn and the vertical wind shear s = ∂u/∂z.
Those parameters were selected because they are known to have a major impact on the CC life-cycle and/or they are expected to affect the relative differences between single-aircraft and formation flight CC. Other parameters that also affect the CC evolution like the initial temperature at flight altitude T * = 217 K are not varied. Those fixed parameters are listed in Table 1. Table 4 gives a list of all parameter variations and combinations. The first block of the section "Variation of the atmospheric scenario" lists all simulations with a final adiabatic cooling ∆T = 2 K where natural cirrus formation by homogeneous nucleation can be neglected and is suppressed in the model. For a moderate wind shear of s = 0.002 s −1 , the updraught speed w syn takes the values 1, 2 and 5 cm s −1 and the initial relative humidity is either RH i,0 * = 110% or 120%. Higher wind shear s = 0.006 s −1 is used in one atmospheric scenario with w syn = 5 cm s −1 and RH i,0 * = 120%.
The second block lists all simulations with a final adiabatic cooling ∆T = 4 K and natural cirrus formation around the CC. Wind shear is fixed at s = 0.002 s −1 , and all four combinations of w syn = 2 or 5 cm s −1 and RH i,0 * = 110% or 120% are tested.
The second table section "Variation of the formation flight scenario" summarises the tests where the FF contrail initialisation is modified. Unterstrasser and Stephan [23] provides data of 3D simulations, where the two aircraft in the formation have a different vertical offset DZ and lateral offset DX. These changes affect the vortex trajectories and with it the early contrail structure. All non-solid red lines in Figure 2 depict the vertical and transverse profiles of simulations with varied DX (45,50,55, and 60 m). Those data are used as input and contrail-cirrus simulations have been performed for one selected atmospheric scenario with RH i,0 * = 110%, w syn = 5 cm s −1 , ∆T = 2 K and Moreover, four different realisations of the same simulation are performed. To do so, the contrail location during the initialisation procedure is horizontally shifted relative to the background turbulence, i.e the location of the black box in Figure 3 changes (see the dotted boxes for the alternative contrail locations). Then, the contrail faces slightly changed turbulence structures which triggers initially small changes in the CC growth and spreading. Those small perturbation may or may not increase over time and the set of four realisations gives a notion of the spread of possible CC evolutions. In the present study, we do not focus on such aspects and only average values of the four realisations are analysed.
Overall, more than 130 CC simulations have been performed to account for the atmospheric variability and the turbulence induced uncertainty. In an FF sensitivity test, further formation flight scenarios are used as input and those simulations are referred to as FF_DX??, as the lateral offset DX of the two aircraft in formation was varied in preceding vortex phase simulations. Various atmospheric scenarios are listed: w01, w02 and w05 refers to an updraught speed w syn = 1, 2 and 5 cm s −1 , s2 and s6 refers to wind shear s = ∂u/∂z = 0.002 or 0.006 s −1 . The initial relative humidity RH i,0 * is either 110% or 120%. Simulations without natural cirrus formation and final adiabatic cooling ∆T = 2 K use the code dt2K and simulations with natural cirrus formation around the CC use the code dT4K. The codes in parentheses or after the single braces starting with "F" refer to the number of the figure (and column) in which the simulation is depicted.

Analysed Quantities
In this final preparatory subsection, we define quantities that will be analysed in the subsequent results section.
The optical thickness τ(x) in a grid column is given by z χdz, where χ(x, z) is the extinction coefficient of the all crystals in a grid box. The total extinction E of a CC is defined as The definition of CC width W OD considers its visibility by a human observer. It counts all CC columns with τ > τ vis = 0.02. The mean optical thickness τ m is the average of τ over all CC columns with τ ≥ τ c = 0.005. The total ice mass I is the integral of the ice mass concentration (IWC) over the contrail cross-section. Analogously, the total IC number N is derived from the IC number concentrations. Hence, N and I have units m −1 and kg m −1 (i.e., per meter of flight path). More information on the definitions can be found in Section 2.3 of Unterstrasser et al. [25].
The climate impact of a specific climate forcing is often measured in terms of a radiation perturbation at the top of the atmosphere (TOA), i.e., radiative forcing or effective radiative forcing [41,42]. Assessing the climate or contrail mitigation potential of formation flight, it would be favourable to evaluate the radiative forcing (RF) of the simulated contrails. However, this quantity cannot directly be derived from the present data as it requires further knowledge on the incident radiation fluxes, which are typically not characterised in LES setups.
Moreover, the contrail RF, which is given in units of W/m 2 , typically quantifies the total radiative effect of all contrails inside a certain area and at a certain moment. It is not an ideal quantity for comparing the radiative effect of individual contrails as it does not take into account changes of the lifetime and its spatial extent. Integrating area-based quantities over the contrail length or horizontal extent and life cycle gives more meaningful quantities for the current approach, in analogy to the energy forcing introduced by Schumann et al. [43,44]. The contrail-intrinsic properties that control its TOA radiative perturbation are the optical thickness τ and the effective radius r eff [45], or equivalently τ and the ice water path IWP, which is the vertial integral of IWC. In analogy to the energy forcing, our proxy metrics for the comparison of the radiative impact of two contrail scenarios are then the time-integrated total extinctionẼ and total ice massĨ.

Exemplary Simulation
First, we present snapshots of an exemplary SA simulation which show the CC cross-sections at initialisation and after 2, 4 and 8 h (Figure 4). The figure displays the extinction coefficient χ. The atmospheric conditions are RH i,0 * = 120%, w syn = 5 cm s −1 , ∆T = 2 K and s = 0.002 s −1 . Here, the simulation time t = 0 h refers to a 5 min old contrail, cf. with Figure 1 of Unterstrasser [33]. The initial contrail is more than 450 m deep and nearly 300 m broad. The contrail is broadest in the lower part, where the Crow instability [46] leads to oscillations of the vortex tube along flight direction (over which is averaged here). Due to the prescribed wind shear, the CC gets tilted and spreads horizontally over time. The largest ICs start to settle and fall out from the tilted CC. A fall streak shapes up, with extinction values that are smaller than in the CC core region, i.e., the thin top layer with much higher χ values. The CC core region features high IC number densities and the fall streak is continuously fed by ICs that keep falling out of the CC core. Below the ISS layer, i.e., z 700 m, ICs shrink until they completely sublimate. Apparently, they disappear prior to reaching the lower domain boundary; we consider those ICs to be lost by sedimentation even though they are actually lost by sublimation. After 2 and 4 h, the width of the fully-grown CC attains values of 9 km and 17 km, respectively. The last snapshot after 8 h shows a decayed CC and the χ values are much smaller. A more detailed discussion of the CC evolution and its response to variations of the atmospheric background can be found in [18,25,26,36,37].

Contrail-Cirrus Evolution
Next, we showcase the variability of the CC evolution across various atmospheric scenarios and also highlight the differences between FF and REF scenarios.
For this, we choose five different quantities to be analysed in more detail. These are the total extinction E, the total ice mass I, the total IC number N , the mean optical thickness τ m and the CC width W OD . Figure 5 depicts only simulations with ∆T = 2 K and no natural cirrus formation. The black vertical bars in the upper row panels indicate the end of the prescribed updraught motion. The first three columns show simulations, where only the updraught speed, and with it the updraught period, is varied. The fourth column uses RH i,0 * = 110% instead of RH i,0 * = 120%; otherwise the parameters are identical to column 3. The fifth column uses a stronger wind shear (s = 0.006 s −1 instead of 0.002 s −1 ); all other simulation parameter are again identical to column 3. We first focus on ubiquitous features of the results, then on differences across the atmospheric scenarios, and finally on the differences between REF and FF scenarios. Several features occur in all five simulations: The total extinction and the ice mass increase as long as the updraught prevails and atmospheric water vapour (WV) is supplied for IC growth. Clearly, after some hours the amount of ice mass is several orders of magnitude higher than the initial contribution from the emitted WV. Soon after the updraught comes to a halt, E and I start to decrease as sedimentation induced losses can no longer be compensated by WV uptake. The peak ice mass often amounts to more than 10 kg per meter of flight path. The IC number and the mean optical thickness decrease monotonically with time, in the first hour often at a much higher rate than later on. In the end, the mean optical thickness is below 0.1.
Next, we discuss the sensitivities to the atmospheric parameters. If w syn is reduced (columns 1 to 3), the peak and final values of E are higher, maximum CC width is larger, more ICs are lost in the first hour, for t > 1h, however, the loss rates are smaller as fewer ICs are lost by sedimentation. If RH i,0 * is reduced (column 3 vs. 4), the lower IC number remains over the total simulation period. Similarly, all other quantities also attain smaller values for smaller RH i,0 * . There are two reasons for this behaviour. First, the amount of supplied WV is smaller by one quarter, as can be derived from the RH i,f * values listed in Table 2  Second, shear-induced horizontal spreading is less effective for an initially shallower CC. Lastly, vertical wind shear is increased (column 3 vs. 5). Note that the strong wind shear simulation is analysed only up to 5.5 h. Checking simulation snapshots of the type shown in Figure 4 reveals that a domain broader than 80 km would be necessary for a longer simulation time. Note that a simple example calculation (∆W = s H ∆t span ) shows that a H = 1 km deep object broadens by around ∆W = 20 km during one hour t span = 1 h. Indeed, we observe that CC width increases much faster in the strong wind shear case. As a consequence, more WV can deposit on the ICs and total ice mass and extinction attain larger values. Moreover, stronger thinning leads to slightly smaller optical thickness values.
Next, we discuss qualitative differences between the REF and the FF cases. The major result is that the Iand E curves of the FF cases lie substantially below the "REF * 2" curves implying strong saturation effects in formation flight scenarios. Often the FF curves are close to the REF curves, suggesting that initial differences in N 0 and CC depth do not matter too much or compensate each other. The saturation effect comes, on first order, from the fact that in a formation two aircraft produce a single contrail that faces the same ambient conditions, in particular experiences the same WV supply, and evolves similarly to an SA CC. Moreover, we find a qualitatively different width evolution berween REF and FF. We first cover the four cases with RH i,0 * = 120% (i.e., all but the fourth column). Within the first two hours the width increase is basically the same for REF and FF. Over the next few hours the REF CCs spread faster, a peak width value is attained after around 4 to 6 h, and width decreases from then on. The FF CCs spread at a lower rate compared to REF from t = 2 h onwards, yet keeps increasing over (nearly) the total simulation time. Often, in the end FF CCs are broader than their REF counterparts.
In order to interpret the width evolution, two aspects have to be considered: First, how far are ICs transported by air motion and sedimentation? This determines the geometric cross-section (defined as the smallest polygon surrounding all CC ICs) and with it the geometric width. However, our width definition is based on visibility, hence the second question is: Are there enough ICs with sufficient ice mass in a specific column to be visible? The REF CCs are initially deeper, hence it seems reasonable that a stronger shear-induced spreading leads to the larger REF width values in the intermediate time period. The mean optical depth decreases over time, and towards the end an increasingly larger number of CC columns becomes sub-visual, hence W OD decreases. In contrast, the geometric width of the FF CCs is smaller and there is less thinning. In combination with a higher IC number, the optical thickness hence remains above the visibility threshold τ vis in most CC columns. This explanation is corroborated by the RH i,0 * = 110% case in column 4. Due to IC loss, the initial SA contrail is shallower than for RH i,0 * = 120% and the contrail depth value is closer to its FF counterpart. Then, we observe a very similar width evolution over the first five hours as the geometric width of FF and REF probably evolve similarly. After five hours, visible CC width keeps increasing only for FF, as the larger optical depth maintains visibility.
Finally, we discuss the NNN simulations, which help to disentangle effects of initial differences in contrail depth and IC number and will reinforce several arguments made above. By construction, NNN and FF simulations start with the same N 0 . Yet, the initial contrail dimensions are equal to REF. Due to the latter, we observe that contrail width of REF and NNN evolves basically identical. Only when the visibility aspect comes into play after several hours, the higher IC numbers and τ values result in a broader visible contrail in the NNN-cases. Moreover, the ice mass evolution of REF and NNN is similar, suggesting similar CC cross-sections and access to ambient supersaturation. In general, light scattering is increased if the same amount of ice mass is distributed over more ICs. For this reason, NNN simulations have higher E values than REF simulations. In most cases, the FF simulation lies in between of NNN and REF. Hence, the E values of the FF simulations are not as high as could be expected from the higher IC number. This implies that the weaker spreading of the FF CC partly offsets this increase in E. This finding helps to refine the above statement that initial differences in N 0 and contrail depth do not matter too much. In fact, they partly cancel out each other and the net effect of both on the saturation extent is indeed of second order.
So far, scenarios with ∆T = 2 K have been discussed. Next, simulations with ∆T = 4 K, where the formation of natural cirrus is allowed, are presented. To shorten the description only E, I and N are depicted in Figure 6. A detailed analysis of how CC and natural cirrus interact on a local scale and how CC evolution is affected by natural cirrus formation is given in Unterstrasser et al. [26] and the reader is referred to this paper for a deeper understanding. Inside the CC, relative humidity is too low to support nucleation of new ICs. ICs form only outside of the CC, once RH i surpasses RH i,crit of nucleation. Natural cirrus surrounds the CC and depletes WV in the vicinity of the CC. This confines the CC growth at its periphery. On the other hand, the sustained ascent continues to make more WV available inside the CC. Hence, in this region IC growth is stronger relative to the 2 K scenarios. These two counteracting effects lead apparently to the result that the Eand I evolutions do not differ much between a 4 K simulation and its respective 2 K counterpart (compare, e.g., first columns in Figures 5 and 6). Below the tilted CC, natural cirrus has been depleting the WV. This weakens the formation of the CC fall streak, as the contrail ice crystal do not grow as fast and sediment more slowly. For this reason, we find the number of ICs to decrease more slowly in the 4 K scenario as fewer ICs are lost by sedimentation.
Comparing the FF CC with REF CC, the findings are very similar to those of the 2 K scenarios. For RH i,0 * = 120% (columns 1 and 2) the FF CC are only slightly stronger in terms of E and I than the Only in the case with RH i,0 * = 110%, the Eand I values of the FF CC are substantially larger. Figure 7 shows again E and I, but now in a normalized version. For this, the FF simulation results are normalized by the "2 * REF" values, e.g, a value of 0.75 means that the "strength" of a contrail produced by an aircraft in a formation is only 75% of that by the same aircraft on a single mission. The added value of the normalized version is manifold. The compact form allows to include the complete set of investigated atmospheric scenarios, simplifies the comparison between FF and REF and nicely summarises conclusions made before. Furthermore, finally, it enables a quantitative comparison.  One apparent finding is that the normalized quantities have larger values (i.e., the saturation effect is less pronounced) for smaller RHi* as the dotted lines lie above the solid lines. In all RH i,0 * = 120% cases (solid lines), the normalized values remain around 0.5 to 0.6 for many hours. In the 4 K cases (right column) it stays like this over the total simulations period. Only in the 2 K cases (left column), the values increase towards the end. This is due to the fact that the FF CC dissolution is slower. However, one has to keep in mind, that the absolute values are already smaller in this stage. Hence, the relative importance of the normalized values on the time-integrated value is somewhat reduced. Now we switch to the RH i,0 * = 110% cases (dotted lines). In the first hour the E norm values drop from 0.85 to 0.65 − 0.70 and then continuously increase over the total simulation period reaching values greater than unity in the end. The I norm values start from 0.5 and increase over time to roughly 1. The CC evolution in absolute numbers depended dominantly s and w syn . Interestingly, the normalized values, which spotlight the differences between REF and FF, are rather unaffected. This implies that REF and FF contrails respond similarly to changes in these atmospheric parameters.
In the first several minutes, the time series of E and I feature small wiggles. Moreover, they suffer from a systematic bias as the FF contrails are a few minutes older than the REF contrail at initialisation. This leads to irregular patterns and overestimates the normalized values E norm and I norm . Those irregularities are unimportant as the absolute values of E and I are still small. For this reason, E norm -and I norm curves in Figure 7 are shown only for t > 30 min.

Time-Integrated Contrail Properties
So far, we assessed time series of contrail properties with the primary goal of explaining the physics behind the saturation effects. In the subsequent and final step of the present study, we analyse single-valued metrics to assess the saturation effect by formation flight. For this, we integrate E and I over time.X where X stands for E or I, and scen for REF, NNN or FF. In order to obtain normalized valuesX,X is divided by the time-integrated "2 * REF" value, i.e., X scen =X scen 2 ·X REF .
The default integration period t int is the total simulation period t sim = 8 h. Assuming that this time period represents a major fraction of the CC life-cycle, the obtained values serve as a measure of the CC radiative impact.
The prescribed ambient humidity evolution favours long contrail lifetimes, in particular for weak, but enduring updraughts. In the present simulations, CC dissolution is mainly triggered by sedimentation. ICs fall into the sub-saturated layer where they sublimate and the ISS layer above becomes dehydrated over time. In reality, ascending air masses may start to subside at some time and, as a result, CCs may start to vanish earlier than in our scenarios. Hence, CC lifetime is limited by sedimentation and/or subsidence as was specifically examined by Bier et al. [47]. To mimic scenarios with smaller contrail lifetimes, we do not perform additional simulations with scenarios that include a downdraught period and which would explicitly simulate the contrail dissolution. Instead, we choose a brute-force approach by analysing the existing simulations and simply shorten the integration period t int in the above formulas. In addition to t sim = 8 h, values of 2 h, 4 h and 6 h are used. The by far largest impact on the saturation ratio has the initial relative humidity. For RH i,0 * = 110% values lie mostly in the interval [0.7, 0.8] and we can see a decrease in saturation for longer contrail lifetimes. Overall, we find that formation flight reduces the time-integrated total extinction by 20% to 50% (corresponding toÊ ∈ [0.5, 0.8]). Compared toÊ, the values of the metricÎ are generally smaller, the parameter RH i,0 * is not as dominant and the sensitivity to t int is larger. Overall, the time-integrated total ice mass is reduced by 30% to 60% (corresponding toÎ ∈ [0.4, 0.7]). In Section 3.2, we mentioned two ways of how a variation of RH i,0 * affects the CC properties: the initial contrail properties and the amount of available ambient WV are changed. In the normalized quantities, the WV effect should, however, cancel out. Hence, the dominant impact of RH i,0 * is due to the fact that the initial contrail properties change with RH i,0 * . This justifies a-posteriori the somewhat detailed presentation of the initial contrail properties given in Section 2.3 and Figure 2. Two aspects of initial contrail properties differ between REF and FF, the geometric cross-section and the total number of ice crystals. The hybrid NNN simulations can help to answer the question which type of initial difference plays a bigger role for the long-lasting differences between FF-and SA CC. NNN  Overall, we conclude that in our scenarios formation flight leads to reductions in time-integrated total extinction and total ice mass by 20% to 50% and 30% to 60%, respectively. This implies a very large mitigation potential of formation flight.

Sensitivity to Formation Flight Geometry
The preceding Section 3.3 demonstrated that the derived saturation effects depend most strongly on the contrail properties after vortex breakup. There, the differences in early contrail properties were triggered by a variation of RH i,0 * . As outlined in Section 2.3, geometric formation parameters like the horizontal offset DX of the two aircraft in formation also have an impact on the early FF contrail properties. So far, however, all presented FF simulations were based on one selected formation geometry with DX = 50 m. Hence, it is pending to check how strongly the derived saturation effects depend on the formation geometry. For this, we select one atmospheric scenario (as detailed in Section 2.4) and initialise the CC simulations with FF contrails with further DX values. Figure 9 shows All four FF scenarios result in qualitatively similar contrail evolutions. Some differences are apparent in the width evolution from t = 6 h onwards and initial differences in N remain over the complete simulation period. Yet, the radiatively important quantities E and I feature only a small spread across the FF scenarios. The panels on the right-hand side depict time series of E norm and I norm in analogy to Figure 7 and confirm the modest to low sensitivity to the FF scenario. Finally, Figure 10 summarizes the above results in the two metrics,Ê andÎ. As in Figure 8, each curve connects four data points, which show values integrated over t int = 2 h, 4 h, 6 h and 8 h, respectively.Î values lie in the range between 0.5 and 0.65. Note that a variation of the atmospheric scenario for fixed RH i,0 * = 110% producedÎ values in the same range. Clearly, the sensitivity to the prescribed integration period t int is stronger than to the FF scenario. In particular for t int = 2 h,Î ∼ 0.52 irrespective of DX. For higher t int , I values exhibit some spread across the FF scenarios. The solid/dotted lines show scenarios with ∆T = 2 K and 4 K, respectively. The set of scenarios is as in Figure 7. Each curve connects four identical symbols and is slightly slanted as the four symbols have a small horizontal offset. For each curve the integration time increases from left to right. The first two columns shows absolute values of total extinction, total ice crystal mass and number, and contrail width (analogous to Figure 5). The right-most column shows normalised values of total extinction and total ice crystal mass (analogous to Figure 7). The grey lines show the reference simulation REF (solid), 2 * REF (dashed) and NNN (dotted). The synoptic scenario is given by w05_dT2K_s2 and RH i,0 * = 110% (see Table 4).
The extinction-based metricÊ reveals some dependence on the FF scenario. In particular, the default DX = 50 m simulation produces largerÊ values than the other DX simulations. The present finding may be a hint that the saturation effects can be larger (i.e., smallerÊ) than what Figure 8 alone suggests as those simulations were solely based on DX = 50 m simulations.

Discussion
In this section we discuss implications and generalisations of our results and also assumptions and limits of the present study.
The saturation effect due to formation flight was evaluated for a multitude of atmospheric scenarios determined by w syn , RH i,0 * , s and ∆T. Those parameters were selected because they are known to have a strong impact on the CC life-cycle and/or they are expected to affect the relative differences between SA and FF CC. For example, the initial contrail depth H shows a characteristic difference between SA and FF scenario. Hence, SA and FF CC may respond differently to a change in vertical wind shear s as CC spreading depends on H and s. In general, we found that the extent of saturation is fairly insensitive to our variations of the atmospheric scenario or the formation geometry. Astonishingly, the initial RH i,0 * was found to be the dominant parameter which deserves a closer investigation. A decrease in RH i,0 * leads to shallower contrails with fewer ice crystals, as most ice crystals in the primary wake sublimate during the early vortex descent. In particular, lowering RH i,0 * from 120% from 110%, as done here, leads to a strong reduction in the initial SA contrail depth see Figure 2 here or Figure 1 in [33]. Contrarily, the FF contrail depth does not change much with RH i,0 * . Thus, the relative differences between the initial SA and FF contrails change substantially by a RH i,0 * change. Other parameters that have not been varied and potentially trigger relative differences between SA and FF cases may have a non-negligible impact on the observed saturation. These are parameters that affect the early contrail evolution and depth such as the aircraft type or the thermal stratification [24]. Our simulation setup, e.g., neglected several processes such as natural cirrus formation by heterogeneous nucleation, contrail radiative heating and aggregation. Even though this may affect the contrail-cirrus evolution, we do not expect that their neglect introduces biases in the comparison of SA and FF scenarios. Similarly, we do not expect that using a 2D-model domain instead of a 3D-model model [37,48] introduces systematic biases.
Moreover, we studied only formations of two identical aircraft and assumed equal emission characteristics for both. In particular, the fuel flow of the follower AC is not reduced, which neglects the reduction of WV emission and the number of generated ICs. Given that a much larger change of ice crystal number between the REF and NNN cases had only modest impact on the amount of saturation, a slight adaptation of fuel flow would result in only marginal changes.
In regions of dense air traffic a "natural" saturation effect occurs when several contrails overlap and form a contrail cluster. In such a case, it is not appropriate to choose an isolated contrail simulation as reference state and the saturation gain due to formation flight might be smaller than the ones we derived. On the other hand, we expect stronger saturation effects in formations with more than two aircraft.
We compare formation and single aircraft contrails at identical atmospheric scenarios. If formations generally fly, e.g., at a lower/higher altitude than the aircraft would do in solo missions, the systematic shift to higher/lower temperatures would affect the contrail properties. Possible biases by such re-routing induced changes of the atmospheric background state are not accounted for here.
The present work is to our knowledge the first study that estimates contrail saturation effect due to formation flight and as such provides a first estimate. It is out of scope to scan the full parameter space of contrail evolution. The implied accuracy of our saturation estimates is acceptable, having in mind, that transferring these numbers into global scale contrail models is done in a simplistic way and thus associated with further uncertainties. The LES-derived saturation values are fed into the non-linear climate response model AirClim [49]. Using realistic aircraft flight trajectories of single mission and formation flight emission scenarios [50] within AirClim provides a first global assessment of the formation flight mitigation potential [50,51].

Summary
This study presents high-resolution simulations of contrail-cirrus that originate from a single aircraft (SA) or a two-aircraft formation (FF). The simulations start with 5 to 7 min old contrails at a time when aircraft wake vortices decayed and the initial contrail data are provided by recent simulation study of early formation flight contrails [23]. The simulated time of 8 h covers large parts of the contrail-cirrus life-cycle.
The analysis focuses on determining the saturation effect that occurs when contrails of multiple aircraft are created in close proximity. Two single-valued metrics are defined in order to compare the contrail-cirrus evolution of SA cases, on the one hand, and FF cases, on the other hand, for a multitude of atmospheric scenarios.
The lifetime-integrated total extinctionẼ and total ice massĨ behind a two-aircraft formation are found to be substantially smaller than for a corresponding case with two separate aircraft and contrails. In our scenarios, formation flight leads to reductions inẼ andĨ by 20% to 50% and 30% to 60%, respectively.
The atmospheric scenarios include variations of the vertical wind shear, initial background relative humidity, the synoptic updraught speed and time period. Even though, variations of those parameters qualitatively change the contrail-cirrus evolution, the relative differences between SA and FF cases are fairly unaffected and the quantified magnitude of reduction seems to be valid for a wide range of situations. Interestingly, the dominant parameter is the initial relative humidity which affects more the early contrail properties than the later transition into contrail-cirrus. A variation of it provokes large relative differences in the early ice crystal number and contrail depth between the SA and FF cases.
Overall, the potential reductions by formation flight are quite substantial. This renders formation flight a promising measure to mitigate the contrail climate impact.
In companion papers of the FORMIC project (FORMation flight: Impact on Climate), the saturation values are fed into the non-linear climate response model AirClim [49]. AirClim also accounts for non-contrail aviation climate effects and was adopted to account for saturation effects occurring during formation flight [51]. Using realistic aircraft flight trajectories of single mission and formation flight emission scenarios [50] within AirClim provides a first global assessment of the formation flight mitigation potential [50,51].