Individual Condensation Trails in Aircraft Trajectory Optimization

Abstract

Contrails are one of the driving contributors to global warming, induced by aviation. The quantification of the impact of contrails on global warming is nontrivial and requires further in-depth investigation. In detail, condensation trails might even change the algebraic sign between a cooling and a warming effect in an order of magnitude, which is comparable to the impact of aviation-emitted carbon dioxides and nitrogen oxides. This implies the necessity to granularly consider the environmental impact of condensation trails in single-trajectory optimization tools. The intent of this study is the elaboration of all significant factors influencing on the net effect of single condensation trails. Possible simplifications will be proposed for a consideration in single-trajectory optimization tools. Finally, the effects of the most important impact factors, such as latitude, time of the year, and time of the day, wind shear, and atmospheric turbulence as well as their consideration in a multi-criteria trajectory optimization tool are exemplified. The results can be used for an arbitrary trajectory optimization tool with environmental optimization intents.

1. Introduction

The term condensation trail (contrail) summarizes a type of human-induced cloud, developed behind aircraft due to condensation of water vapor emissions and ambient humidity around exhausted soot particles and atmospheric condensation nuclei [1] in a cold ambient atmosphere, satisfying the Schmidt–Appleman criterion [2,3]. In an ice supersaturated environment, these artificial ice clouds form into long-living cirrus clouds called “Cirrus homogenitus” as defined by the World Meteorological Organization [4].

In the Earth–atmosphere energy budget, contrails act like a barrier [5,6,7]. They scatter incoming shortwave solar radiation back to the sky (resulting in a cooling effect) and they absorb and emit the outgoing longwave terrestrial radiation back to the Earth’s surface (yielding a warming effect in the lower layer of the atmosphere) [6,8,9,10]. The dominating effect may be defined as radiative forcing $RF$, as the change in net downward radiative flux at the tropopause (considering the instantaneous response of the stratosphere) [6]. The contribution of contrails to global warming cannot yet be determined with high precision and depends on flight performance, weather conditions, and time. Latest combinations of several approaches to model a global impact of contrails on global warming summarize a warming net effect of $R{F}_{\mathrm{Contrail}}=0.05$ W m${}^{-2}$ with accuracy to within $-0.02$ and $+0.15$ W m${}^{-2}$ for 2010 [6]. Single studies result in a precision of the environmental impact of contrails including negative (cooling) net effects (e.g., $R{F}_{\mathrm{Contrail}}=-0.007\mathrm{to}+0.02$ W m${}^{-2}$ for 2005 [9]). Using global climate models and historic air traffic data, reliable estimations of the global contrail radiative forcing for the year 2000 of $R{F}_{\mathrm{Contrail}}=0.03(-0.01\mathrm{to}+0.08)$ W m${}^{-2}$ [9] have been improved for 2010 to $R{F}_{\mathrm{Contrail}}=0.02(-0.01\mathrm{to}+0.03)$ W m${}^{-2}$ [6] considering an increased traffic distance by 22% between 2005 and 2010. For 2002, Burkhardt and Kärcher [10] estimated $R{F}_{\mathrm{Contrail}}=0.03$ W m${}^{-2}$ of contrails and contrail cirrus clouds within a global climate model.

Compared to the summed effect of one year’s aviation-emitted carbon dioxide ${\mathrm{CO}}_2$ and nitrogen oxides ${\mathrm{NO}}_{\mathrm{x}}$ which is $R{F}_{{\mathrm{CO}}_2,{\mathrm{NO}}_{\mathrm{x}}}=0.04$ W m${}^{-2}$ [6], the impact of contrails on global warming is no longer negligible [6,11], as it might exceed those effects of ${\mathrm{CO}}_2$ and ${\mathrm{NO}}_{\mathrm{x}}$, even though contrails are only formed during $10\%$ of the flights, on a global average [12]. Thus, the need for contrails to be considered in trajectory optimization becomes indisputable.

The impact of single condensation trails on trajectory optimization has been analyzed by Gounou et al. [13] and Forster et al. [14], focusing on the importance of large solar zenith angles during sunset and sunrise. In an application of a Monte Carlo code for photon transport, Forster et al. [14] already considered effects like multiple scattering, but in a coarse spatial grid and ignoring the impact of flight performance on the optical properties of the contrail. Detailed studies on all significant impact factors have been elaborated by Rosenow [11]. In the current study, the results have been simplified and harnessed for trajectory optimization. All other research interests known to the authors concentrated on the effect of contrails on a global scale.

Using satellite data, the adjusted forcing $AF$ as the imbalance of the Earth–atmosphere energy system after stratospheric temperatures have been adjusted to regain a radiative equilibrium in the stratosphere (assuming zero further radiative heating rates) can be calculated [15,16], considering a completed transition of contrails into cirrus. Herewith, the diurnal cycle of contrails and cirrus and differences in regions with low and high air traffic demand can be distinguished [16]. Using satellite data of 2006, $A{F}_{\mathrm{Contrail}}=0.045\mathrm{to}0.075$ W m${}^{-2}$ has been quantified for contrails and contrail-induced cirrus. By applying a combination of modeled and satellite data-based estimates with tolerances in spreading rate, contrails optical depth, ice-particle shape, and radiative transfer [17] and by accounting for the ongoing increase in air traffic, contrails, and contrail-induced cirrus $A{F}_{\mathrm{Contrail}}$ for the year 2010 of $A{F}_{\mathrm{Contrail}}=0.05\left(0.02\mathrm{to}0.15\right)$ W m${}^{-2}$ is widely accepted [6]. Bock and Burkhardt [18] predicted a global future contrail radiative forcing of $R{F}_{\mathrm{Contrail}}=0.159$ W m${}^{-2}$ for 2050 assuming $R{F}_{\mathrm{Contrail}}=0.049$ W m${}^{-2}$ in 2006 and considering a shift in cruising altitudes, an increase in air traffic distance until 2050 by a factor of 4 compared to 2006, an increase of alternative fuels (with reduced soot emission) and expected changes in propulsion efficiency in 2050. They determined the distribution of the contrail formation taking into account the global distribution of the expected air traffic distance. However, the influence of those contrails on global warming was again determined using global averages. Chen and Gettelman [19] modeled the impact of existing contrails on contrail cirrus ice-crystal shape and size. They predicted an increase in contrail cirrus radiative forcing by a factor of 7 for 2050 (i.e., $R{F}_{\mathrm{Contrail}}=0.087$ W m${}^{-2}$) compared to 2006 assuming, that air traffic movements will grow on average by a factor of 4 until 2050. In their model, Chen and Gettleman considered a region-specific contrail radiative forcing and a local increase in air traffic. Avila and Sherry [20] estimated the radiative forcing of individual contrails, applying a model developed by Schumann et al. [21], assuming a constant optical depth. Here, contrails are roughly parameterized in terms of optical depth, width, and particle diameter. For each contrail class, an effective particle radius is approximated according to Schumann and each effective radius corresponds to specific optical properties. The large advantage of the study by Avila and Sherry is the consideration of the solar zenith angle.

A consideration of contrails in trajectory optimization tools, if any, has only been found by taking into account a constant value for trajectories through ice-supersaturated regions [22,23,24,25,26]. This consequent avoidance of contrail formation does not, however, lead to a holistically optimized trajectory [27]. Specifically, considering single contrails with a cooling effect on the Earth–atmosphere system is completely misunderstood in those approaches. The rest of the paper is organized as follows: In Section 2, the methodology of each sub-model is presented, whilst in Section 3 the implementation of all sub-models in the trajectory optimization is demonstrated. Furthermore, a validation of the model and a brief discussion of the findings is shown. In Section 4, the findings have been generalized for an easy implementation in a fast-time trajectory optimization tool.

2. Individual Contrails in Trajectory Optimization

The assessment of single contrails in trajectory optimization may be performed in six steps. The procedure is shown in Figure 1.

2.1. Atmosphere GFS Weather Data

First, detailed weather information for the environment around the aspired route is required. Modeled and gridded global forecast data from the Global Forecast System (GFS), provided every six hours in Grib2 format by the National Oceanic and Atmospheric Administration (NOAA) with a spatial resolution of $0.25$ degrees at 18 pressure levels is a good compromise between computational effort and accuracy. For trajectory optimization and the estimation of the conditions of contrail formation, vertical profiles of temperature T [K], pressure p [Pa], density $\rho$ [$\mathrm{kg}{\mathrm{m}}^{-3}$], horizontal wind components u [$\mathrm{m}{\mathrm{s}}^{-1}$] and v [$\mathrm{m}{\mathrm{s}}^{-1}$], vertical velocity w [$\mathrm{Pa}{\mathrm{s}}^{-1}$], and the relative humidity $rH$ [-] along the whole route are taken into account.

2.2. Flight Performance

Together with typical input variables for trajectory optimization, such as city pair, aircraft type, engine type, payload, optimization function (i.e., minimum fuel burn, minimum time of flight, minimum contrail impact, or multi-criteria optimization), a trajectory optimization model with implemented key performance assessment can be used for the calculation of the optimum vertical and lateral path. Here, we use the validated simulation environment TOMATO [28,29,30] which includes the aircraft performance model COALA [31,32] for vertical optimization and for the quantification of the emissions. In TOMATO, the trajectory is optimized iteratively by assessing each interim solution regarding several key performance indicators (KPI) including contrails [33]. For comparability, each KPI is transformed into monetary values using the European Emission Scheme. The weighting of the contrail impact in the trajectory optimization is adjusted with prices between €10 and €65 per tonne CO${}_2$ equivalent emission. According to the Institute for Advanced Sustainability Studies, a price of €65 per tonne CO${}_2$ equivalent emission is required to restrict global warming to a limit of $dT=2$ °C [34].

For the evaluation of the contrail, two conditions on contrail formation have to be fulfilled: Firstly, the Schmidt–Appleman criterion [2,3] determines the threshold temperature $T_{\mathrm{LC}}$ of the atmosphere which is required to overcome the exhaust heat. $T_{\mathrm{LC}}$ depends on the aircraft true air speed $v_{\mathrm{TAS}}$ [$\mathrm{m}{\mathrm{s}}^{-1}$], thrust F [N], fuel flow ${\dot{m}}_f$ [$\mathrm{kg}{\mathrm{s}}^{-1}$], specific combustion heat of kerosene $Q=43\mathrm{MJ}{\mathrm{kg}}^{-1}$, and on the emission index of water vapor $E{I}_{\mathrm{water}}=1.24{\mathrm{kg}}^{-1}$ [11]. The mixing of the exhaust emissions with the ambient atmosphere follows a straight line in an enthalpy-mass diagram with the slope G

$$G=\frac{Q(1-\eta )}{E{I}_{\mathrm{water}}},$$

(1)

where $\eta$ defines the overall propulsion efficiency as the fraction between resulting energy (thrust F and true air speed $v_{\mathrm{TAS}}$) and required energy (specific combustion heat Q and fuel flow ${\dot{m}}_f$)

$$\eta =\frac{F{v}_{\mathrm{TAS}}}{Q{\dot{m}}_f}.$$

(2)

In case of the mixing line of the exhaust emissions with the ambient atmosphere touching the saturation water vapor pressure, contrail formation takes place (see [2,3,11] for more details). From this, it follows that the larger the G the higher the threshold temperature for contrail formation and the higher the probability of contrail formation. The second criteria for contrail formation is an ice-supersaturation of the ambient atmosphere in order to survive the vortex regime of the wake vortices, wherein the exhaust emissions are captured [35,36]. Therewith, sufficient humidity for a constant ice-particle growth is provided. Each criteria, the Schmidt–Appleman criterion and the ice-supersaturation, is examined for each second of the flight performance calculation.

2.3. Contrail Life Cycle

In the case of both criteria, the Schmidt–Appleman criterion and the ice-supersaturation, being satisfied, T, p, $\rho$, u, v, w, $rH$, longitude $lon$ [°], latitude $lat$ [°], altitude z [m], ${\dot{m}}_f$ [$\mathrm{kg}{\mathrm{s}}^{-1}$], and time t [s] are provided by COALA to the “Contrail Life Cycle” model [11,37] for each time step, a contrail is induced. For the initial dimensions, a wake vortex model [38] is applied and calibrated for each aircraft type, implemented in COALA. The initial decay of the contrail in the vortex regime strongly depends on atmospheric turbulence (see Figure 2, right), which is approximated by the energy dissipation rate $\epsilon$ [${\mathrm{m}}^2{\mathrm{s}}^{-3}$] as conversion of kinetic energy due to molecular friction per unit mass and per time into thermal energy [39]. $\epsilon$ is calculated assuming a lognormal distribution of turbulence in the lower troposphere and upper stratosphere and a linear correlation between a logarithmic diagnostic turbulence value, such as the vertical velocity w [Pa s${}^{-1}$], provided by the GFS [40,41]. The vertical velocity w is defined as the Lagrangian rate of change of pressure with time [42]

$$w=\frac{dp}{dt}.$$

(3)

Furthermore, wind shear $sh$ [s${}^{-1}$] (Figure 2, left) as difference in wind velocity u and v [$\mathrm{m}{\mathrm{s}}^{-1}$] between two altitudes $\Delta z$ [m] strongly influences the 2D sheared Gaussian plume model for contrail dispersion [11,37]. $\Delta z$ is called shear layer and depends on the maximum differences in wind velocity between two altitudes. $sh$ is also calculated from the provided weather data. The sheared diffusivity $D_s$ [m${}^2$ s${}^{-1}$] is assumed to be in the range of the square root of the vertical $D_v$ [m${}^2$ s${}^{-1}$] and horizontal diffusivity $D_h$ [m${}^2$ s${}^{-1}$]: $D_s\approx \sqrt{D_v{D}_h}$ [43], but $D_s\le \sqrt{D_v{D}_h}$ [11]. Assuming a soot emission index of $E{I}_{\mathrm{soot}}=0.04\mathrm{g}{\mathrm{kg}}^{-1}$ kerosene [44] and a proportional share of ice-particles in the contrail [44], all variables impacting the optical properties of the contrail can be provided to the “Contrail Life Cycle” model. In fact, these are the horizontal ${\widehat{\sigma }}_h\left(t\right)$ [${\mathrm{m}}^2$], vertical ${\widehat{\sigma }}_v\left(t\right)$ [${\mathrm{m}}^2$], and sheared ${\widehat{\sigma }}_s\left(t\right)$ [${\mathrm{m}}^2$] components of the contrail diffusivity variance $\widehat{\sigma }\left(t\right)$ [${\mathrm{m}}^2$], the contrail cross section $CCS\left(t\right)$ [m${}^2$], the ice-water content $IWC\left(t\right)$ [$\mathrm{kg}{\mathrm{m}}^{-3}$] as total amount of ice mass per volume contrail, the number of ice-particles $n_p$ (hereafter called ice-particle number density) [${\mathrm{m}}^{-3}$] and the projected particle area $A_p$ [m${}^2$]. The impact of $sh$, $\epsilon$, and vertical wind speed $v_z$ [$\mathrm{m}{\mathrm{s}}^{-1}$] on the particle radius and the contrail lifetime is shown in Figure 2 and Figure 3, respectively. The vertical wind speed $v_z$ can be calculated from the mean divergence of the horizontal wind speed, averaged along the vertical axis between ground and flight level [45]. Furthermore, $v_z$ can be calculated from vertical velocity w using

$$v_z=\frac{w}{\rho g}$$

(4)

where $\rho$ [kg m${}^{-3}$] denotes the ambient air density and g [m s${}^{-2}$] the gravitational acceleration. Thereby, negative values indicate upward wind speeds. On average, values of $\overline{v_z}=-0.005$ m s${}^{-1}$ are calculated. The vertical wind velocity is two orders of magnitudes smaller than the horizontal wind velocity and a vertical upwind velocity of $v_z=-0.005$ m s${}^{-1}$ is realistic in stable stratification [46]. Assuming this vertical upward speed and a sedimentation of the ice-particles following Stoke’s law, the contrail sediments as soon as the ice-particle radius exceeds values of $r_p\approx 6\mathsf{\mu }\mathrm{m}$ [11].

Note, the available ice-water content is distributed equally to the available number of ice-particles per contrail volume. During cruise, aircraft fly with a speed between maximum Mach number (upper boundary) and a maximum specific range (lower boundary). Within this range, fuel flow decreases with decreasing speed. Since the number of emitted soot particles is assumed to be proportional to the fuel flow, a decreasing fuel flow (as a consequence of a decreasing speed) causes a lower number of soot particles (as primary condensation nuclei). The water vapor in the atmosphere and from the emissions thus sublimates to a lower number of soot particles. As a result, fewer but larger ice-particles are formed. However, the optical properties of contrails depend more strongly on ice-particle size than on ice-particle number density [11]. From this follows that flying slowly through ice-supersaturated regions causes more radiative effective contrails. Furthermore, according to the Schmidt–Appleman criterion [2,3], the higher the overall engine efficiency, the lower the exhaust gas temperature, the higher the critical atmospheric temperature ($T_{\mathrm{LC}}$) for contrail formation [47]. The critical atmospheric temperature for contrail formation $T_{\mathrm{LC}}$ positively correlates with thrust, true air speed $v_{\mathrm{TAS}}$, and fuel flow ${\dot{m}}_f$. The higher the $T_{\mathrm{LC}}$, the higher the probability of contrail formation in an ice- supersaturated region. In general, optimized air speeds for highly efficient conditions of combustion (i.e., minimum fuel flow) refer to low values of $TAS$. Both facts support the thesis that cold, ice-supersaturated regions should preferably be flown through with high speeds.

2.4. Contrail Optical Properties

The fourth step is the “Contrail Optical Properties” model. Depending on the geometrical and microphysical properties of the contrail, the radiative extinction due to scattering, absorption, and emission within the contrail is calculated by running a Monte Carlo Simulation to consider multiple scattering events which are likely, especially for large solar zenith angles $\theta$ [rad] [11]. Therefore, Beer’s law

$$\frac{I_{\lambda }\left(s_1\right)}{I_{\lambda }\left(s_2\right)}=exp\left[-{\int }_{s_1}^{s_2}-Q_e{A}_p{n}_p\left(s\right)\mathrm{d}s\right]$$

(5)

is used, where $I_{\lambda }\left(s_2\right)$ and $I_{\lambda }\left(s_1\right)$ denote the original and the extinguished wavelength specific radiation of solar intensities [$\mathrm{W}{\mathrm{m}}^{-2}{\mathrm{sr}}^{-1}$] and of terrestrial irradiances [$\mathrm{W}{\mathrm{m}}^{-2}$] (see Figure 4). $Q_e\left(s\right)$ denotes the extinction efficiency [-] and depends on wavelength, particle size, and shape [11,48] and $A_p$ denotes the projected particle area [m${}^2$]. $Q_e\left(s\right)$ and $A_p$ (s) are not constant within the contrail and depend on the position s. In the “Contrail Optical Properties” model, $\frac{I_{\lambda }\left(s_1\right)}{I_{\lambda }\left(s_2\right)}$ Equation (5) is interpreted as a number ratio of extinguished photons, regardless of the amount of radiation irradiating the contrail (see Figure 4) [11,49].

The extinction of photons is calculated for each direction in space individually, considering ice-particle shape-dependent and wavelength-dependent probability of an extinguishing event. The latter is described by the absorption and scattering efficiency $Q_a\left(s\right)$, $Q_s\left(s\right)$. $Q_s\left(s\right)+Q_a\left(s\right)=Q_e\left(s\right)$ are parameterized by Wyser et al. [48] and Yang et al. [50] as functions of wavelength, ice-particle size, shape, and density, which in turn are provided by the Gaussian plume model (see Figure 4).

This “Contrail Optical Properties” model provides weighted number ratios $S_i(\lambda ,\mathrm{d}\omega )$ of extinguished photons per meter contrail, per wavelength and per time step of the contrail life cycle (Figure 4). For each direction of incoming photons, $S_i(\lambda ,\mathrm{d}\omega )$ [m] is calculated by

$$S_i(\lambda ,t,\mathrm{d}\omega )=\frac{N_{\mathrm{out}}}{N_{\mathrm{in}}}{w}_{\mathrm{in}}sin\alpha ,$$

(6)

where $N_{\mathrm{out}}$ and $N_{\mathrm{in}}={10}^7$ denote the number of outgoing and incoming photons, respectively. $w_{\mathrm{in}}=6{\widehat{\sigma }}_h\left(t\right)$ denotes the irradiated width of the contrail, $\alpha$ defines the angle between the length axis of the contrail and the incoming photons (see Figure 5). This extinction strongly depends on the direction of irradiation (i.e., $\alpha$). The longer the travel distance of photons through the contrail (i.e., the larger $\alpha$), the higher the probability that an extinguishing event takes place. With increasing travel distance and with increasing $\alpha$, the number of out-scattered photons increases, although a dominating forward scattering is expected [48,50,51]. From this follows that during horizontal photon transport during sunrise and sunset, more photons will be scattered in the opposite hemisphere compared to a vertical photon transport at noon. However, the contrail radiative extinction further depends on the power with which the contrail is irradiated, which will be maximum at noon and minimum at night. This power is calculated separately in the next step.

2.5. Atmospheric Radiative Transfer

The fifth step of the approach provides wavelength and angular specific solar intensities and terrestrial irradiances at the position of the contrail. This step is necessary, because the radiative extinction due to the contrail does not only depend on the optical properties of the contrail itself but also on the amount of radiation irradiating the contrail. The powers of solar intensities and terrestrial irradiances depend on wavelength, longitude, latitude, altitude, the presence of clouds, time of the day, and time of the year. Most of these input variables are provided by the flight performance model (see Figure 1). For the calculation of this “Atmospheric Radiative Transfer”, the radiative transfer software package libRadtran [52] is used. Due to different properties of solar and terrestrial extinction in the atmosphere, different radiative transfer solvers are used.

2.5.1. Terrestrial Radiative Transfer

In the terrestrial wavelength spectrum ($3\le \lambda \le 100\mathsf{\mu }\mathrm{m}$), absorption by atmospheric molecules strongly depends on wavelength and varies between neighboring wavelengths. These narrow absorption bands require a high spectral resolution in the radiative transfer calculation and therefore a high computational effort. However, a weak angular dependence (described by zenith angle $\theta$ and azimuthal angle $\varphi$) of radiation due to a missing part of direct irradiance is expected. The two stream approximation (TSA) takes advantage of the weak angular dependence and reduces the computational effort [53]. Here, all shares of radiation coming from one hemisphere are azimuthally averaged over the half-space and are treated as a single irradiance F [$\mathrm{W}{\mathrm{m}}^{-2}$] without specific information conserning the angular direction (see Figure 6, left). Due to this average, two irradiances at any altitude remain: terrestrial irradiances coming from the lower hemisphere

$$F_{\mathrm{up}}(\lambda ,t,\theta ={}^1{/}_2\pi ...{}^3{/}_2\pi )$$

(7)

and terrestrial irradiances coming from the upper hemisphere

$$F_{\mathrm{down}}(\lambda ,t,\theta ={}^3{/}_2\pi ...{}^1{/}_2\pi )$$

(8)

(see Figure 6, left) [11].

2.5.2. Solar Radiative Transfer

In the solar wavelength spectrum ($0\le \lambda \le 4\mathsf{\mu }\mathrm{m}$, see Figure 7) a TSA is out of the question because of a large influence of the direct beam (in the direction of the position of the sun) on the radiation field. This influence causes a strong angular radiative dependence which cannot be described by a TSA. The radiative transfer solver DISORT (DIScrete Ordinate Radiative Transfer solver) is used for the angular-dependent calculation of direct solar intensities [$\mathrm{mW}{\mathrm{sr}}^{-1}{\mathrm{m}}^{-2}{\mathrm{nm}}^{-1}$].

$$I_{\mathrm{dir}}(\lambda ,t,\mathrm{lon},\mathrm{lat},\mathsf{\Omega })$$

(9)

Here, the direction of the direct beam, $\mathsf{\Omega }$ [sr], is described by the solar zenith angle and the solar azimuthal angle. Diffuse solar intensities (Equation (10)) [$\mathrm{mW}{\mathrm{sr}}^{-1}{\mathrm{m}}^{-2}{\mathrm{nm}}^{-1}$]

$$I_{\mathrm{diff}}(\lambda ,t,\mathrm{lon},\mathrm{lat},\theta ,\varphi )$$

(10)

depending on longitude, latitude, altitude, time of the day, and time of the year [54] are also provided by DISORT. DISORT is the most used, recommended and most updated solver for angular-depending radiative transfer in the shortwave spectrum [11,55]. Diffuse solar intensities $I_{\mathrm{diff}}$ are calculated with an angular discretization of $\mathrm{d}\theta =\mathrm{d}\varphi =2$ and the direct beam $I_{\mathrm{dir}}$ is added to the solid angle $\mathrm{d}\omega$ [sr]

$$\mathrm{d}\omega =sin\theta \mathrm{d}\theta \mathrm{d}\varphi ,$$

(11)

where $\theta$ and $\varphi$, are the solar zenith and azimuthal angle, respectively (see Figure 6, right) [11]. Figure 7 clearly indicates the different contributions of wavelength-specific irradiances on the radiation budget of the earth–atmosphere system. Hence, a consideration of solar irradiances around the maximum at $\lambda =0.55$ (i.e., $0.2<\lambda <1\mathsf{\mu }$m) would be sufficient. Terrestrial irradiances should be considered between $5<\lambda <22\mathsf{\mu }$m. Additionally, Figure 7 proves the significant impact of atmospheric absorption by comparing the modeled irradiances with the theoretical values, calculated with Planck’s function assuming mean temperatures of 5750 K of the sun and 288 K of the Earth’s surface.

2.6. Contrail Radiative Forcing

The radiative quantities of the “Atmospheric Radiative Transfer” model are multiplied with the extinguished number ratios, provided by the “Contrail Optical Properties” model to estimate the wavelength-specific and direction-specific extinction of radiation due to the contrail. Because the contrail radiative forcing is defined as imbalance of the radiation budget at a horizontal layer, we distinguish between sources and drains of radiation in the upper and in the lower hemisphere of the contrail (see Figure 8).

To calculate the radiative forcing per unit length of the contrail, the power $P_i(\lambda ,t,\mathrm{lon},\mathrm{lat},\mathrm{d}\omega )$ [$\mathrm{W}{\mathrm{m}}^{-1}{\mathrm{nm}}^{-1}$] of the extinguished photons once irradiated on a unit length contrail have to be considered and balanced. Therefore, the solar intensities in Equation (10) and the terrestrial irradiances in Equations (7) and (8) coming from a particular solid angle $\mathrm{d}\omega$ [sr], (calculated by the “Atmospheric Radiative Transfer” model) have to be weighted by the weighted number ratios of extinguished photons $S_i$ [m] (Equation (6)) (from the “Contrail Optical Properties” model) and the corresponding solid angle $\mathrm{d}{\omega }_i$. Therefore, only the backscattered and absorbed powers are of interest and the extinguished powers can be summarized to the following components (see Figure 8):

• $P_{\uparrow \downarrow a}$: number of absorbed photons coming from below;
• $P_{\uparrow b}$: number of photons coming from below, scattered into the lower hemisphere;
• $P_{\downarrow b}$: number of photons coming from above, scattered into the upper hemisphere.

For radiation coming from a particular solid angle $\mathrm{d}{\omega }_i$ (Equation (11)) and getting scattered into the same hemisphere (i.e., backscattered), the extinguished power $P_b(\lambda ,t,\mathrm{lon},\mathrm{lat},\mathrm{d}{\omega }_i)$ [$\mathrm{W}{\mathrm{m}}^{-1}{\mathrm{nm}}^{-1}$] is

$$\begin{array}{ccc}\hfill P_{\downarrow b}(\lambda ,t,\mathrm{d}{\omega }_i)& =& I(\lambda ,t,\mathrm{d}{\omega }_i)·S_b(\lambda ,t,\mathrm{d}{\omega }_i)·\mathrm{d}{\omega }_i\hfill \\ \hfill P_{\uparrow b}(\lambda ,t,\mathrm{d}{\omega }_i)& =& I(\lambda ,t,\mathrm{d}{\omega }_i)·S_b(\lambda ,t,\mathrm{d}{\omega }_i)·\mathrm{d}{\omega }_i\hfill \\ \hfill P_{\downarrow \uparrow a}(\lambda ,t,\mathrm{d}{\omega }_i)& =& I(\lambda ,t,\mathrm{d}{\omega }_i)·S_a(\lambda ,t,\mathrm{d}{\omega }_i)·\mathrm{d}{\omega }_i,\hfill \end{array}$$

where $S_b(\lambda ,t,\mathrm{d}{\omega }_i)$ denotes the weighted number ratio of backscattered photons coming from a particular solid angle $\mathrm{d}{\omega }_i$ and getting scattered into the same hemisphere and $I(\lambda ,t,\mathrm{d}{\omega }_i)$ is the wavelength-specific solar intensity from the particular solid angle $\mathrm{d}{\omega }_i$. In the same way, absorbed solar powers $P_{\downarrow \uparrow a}(\lambda ,t,\mathrm{lon},\mathrm{lat},\mathrm{d}{\omega }_i$) are calculated for each wavelength $\lambda$ and each solid angle $\mathrm{d}{\omega }_i$.

Terrestrial irradiances $F_{\mathrm{up}}$ and $F_{\mathrm{down}}$ [$\mathrm{W}{\mathrm{m}}^{-2}{\mathrm{nm}}^{-1}$], calculated with the two stream approximation, are hemispherically averaged irradiances. To estimate the extinguished power in the terrestrial wavelength spectrum, $F_{\mathrm{up}}$ and $F_{\mathrm{down}}$ have to be weighted by the number ratio of extinguished photons coming from the upper hemisphere $S_{\downarrow b}$ and the lower hemisphere $S_{\uparrow b}$ with both getting scattered into the same hemisphere that they are coming from and the absorbed photons coming from all directions $S_{\downarrow \uparrow a}$. For example:

$$\begin{array}{cc}\hfill P_{\downarrow b}(\lambda ,t)& =F_{\mathrm{down}}{S}_{\downarrow b}(\lambda ,t)\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill P_{\uparrow b}(\lambda ,t)& =F_{\mathrm{up}}{S}_{\uparrow b}(\lambda ,t)\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill P_{\downarrow \uparrow a}(\lambda ,t)& =F_{\mathrm{down}}(\lambda ,t)S_{\downarrow a}(\lambda ,t)\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill & +F_{\mathrm{up}}(\lambda ,t)S_{\uparrow a}(\lambda ,t).\hfill \end{array}$$

(15)

The backscattered powers $P_{b,}(\lambda ,t,\mathrm{lon},\mathrm{lat},\mathrm{d}\omega )$ are categorized into two classes: Firstly, backscattered powers from the upper hemisphere $P_{\downarrow b}(\theta =\frac{3}{2}\pi ...\frac{1}{2}\pi$) (resulting in a cooling effect) and secondly, backscattered powers from the lower hemisphere $P_{\uparrow b}(\theta =\frac{1}{2}\pi ...\frac{3}{2}\pi$) (with a heating effect, see Figure 8). Absorbed powers $P_{\downarrow \uparrow a}(\lambda ,t,\mathrm{lon},\mathrm{lat},\mathrm{d}\omega$) coming from both hemispheres always contribute to a contrail heating. Those components contribute to the radiative forcing $R{F}_{\mathrm{Contrail}}$

$$R{F}_{\mathrm{Contrail},\mathrm{m}}=P_{\downarrow \uparrow a}-P_{\uparrow b}+P_{\downarrow b}.$$

(16)

After integrating over all significant wavelengths, time steps of the contrail life cycle and over the contrail length, each with different optical properties $S_i(\lambda ,t)$, the total radiative forcing $R{F}_{\mathrm{Contrail},\mathrm{m}}$ [W m${}^{-1}$] (Equation (16)) per meter of an individual contrail can be estimated. Multiplying Equation (16) with the length $L_{\mathrm{Contrail}}$ of the contrail with similar optical properties offers the extinguished power [W] due to the contrail. Usually, the radiative forcing is related to the Earth’s surface, which is $A_{\mathrm{Earth}}=5.1\times {10}^{14}{\mathrm{m}}^2$. Hence, the estimated extinguished power is divided by $A_{\mathrm{Earth}}$.

$$R{F}_{\mathrm{Contrail}}=\frac{R{F}_{\mathrm{Contrail},\mathrm{m}}{L}_{\mathrm{Contrail}}}{A_{\mathrm{Earth}}}.$$

(17)

2.7. Weighting of Contrail Costs in Trajectory Optimization

The integration of $R{F}_{\mathrm{Contrail}}$ over the whole life cycle, over all significant wavelengths and the division by $A_{\mathrm{Earth}}$ offers a comparison with the impact of a reference emission (e.g., ${\mathrm{CO}}_2$) over the whole flight.

Usually, a certain time horizon H (e.g., $H=100$ years) is considered in these interpretations. In this case, the ratio of $R{F}_{\mathrm{Contrail}}$ over $R{F}_{{\mathrm{CO}}_2}$ (known as the global warming potential, $GWP$)

$$GWP=\frac{\underset{0}{\overset{\mathrm{H}}{\int }}R{F}_{\mathrm{Contrail}}\left(t\right)\mathrm{d}t}{\underset{0}{\overset{H}{\int }}R{F}_{{\mathrm{CO}}_2}\left(t\right)\mathrm{d}t},$$

(18)

is used to calculate ${\mathrm{CO}}_2$ equivalent emissions, which are transferred into monetary values by applying political instruments, such as the emission trading scheme (ETS). For instance, for $H=100$ years, the radiative forcing of the total amount of ${\mathrm{CO}}_2$ in the atmosphere is [56]

$$\underset{0}{\overset{H=100}{\int }}R{F}_{{\mathrm{CO}}_2}\left(t\right)\mathrm{d}t={\alpha }_{{\mathrm{CO}}_2}ln\left(\frac{\mathrm{C}}{{\mathrm{C}}_0}\right)=1.94\mathrm{W}{\mathrm{m}}^{-2},$$

(19)

where ${\alpha }_{{\mathrm{CO}}_2}=5.35$ denotes a constant [56], $\mathrm{C}=399.39\mathrm{ppm}$ is the actual concentration of ${\mathrm{CO}}_2$ [57], and ${\mathrm{C}}_0=278\mathrm{ppm}$ is the pre-industrial concentration of ${\mathrm{CO}}_2$ in 1959 [57].

Equation (19) considers the total sum of ${\mathrm{CO}}_2$ emissions in 100 years and is not restricted to the aviation-transport sector. The radiative impact of ${\mathrm{CO}}_2$ strongly depends on the altitude of emission. In high altitudes (above the tropopause with low vertical exchange), ${\mathrm{CO}}_2$ exists for 400 years. In low altitudes (below the tropopause) precipitation induces washing-off effects of ${\mathrm{CO}}_2$. Those large differences in the residence time complicate the use of Equation (19) for a comparison of aviation-induced radiative forcing due to contrails and ${\mathrm{CO}}_2$ emissions.

However, from the literature review in the introduction of this paper, we know the amount of ${\mathrm{CO}}_2$ emissions caused by aviation in 2005

$$m_{{\mathrm{CO}}_2,2005}=733\times {10}^6\mathrm{t}{\mathrm{a}}^{-1},$$

(20)

Ref. [7] and the radiative forcing of those emissions

$$R{F}_{{\mathrm{CO}}_2,2005}=0.028\mathrm{W}{\mathrm{m}}^{-2}.$$

(21)

Hence, the impact of each tonne ${\mathrm{CO}}_2$ emitted by aviation in 2005 was

$$R{F}_{{\mathrm{CO}}_2}=3.8\times {10}^{-11}\mathrm{W}{\mathrm{m}}^{-2}{\mathrm{t}}_{{\mathrm{CO}}_2}^{-1}.$$

(22)

Therewith, the radiative forcing of an individual contrail, calculated in Equation (17), can be transformed into tonnes of ${\mathrm{CO}}_2$ equivalent emissions, which can be used as external costs in a multi-criteria trajectory optimization.

$$m_{{\mathrm{CO}}_2\mathrm{eq}}=\frac{R{F}_{\mathrm{Contrail}}}{3.8\times {10}^{-11}\mathrm{W}{\mathrm{m}}^{-2}{\mathrm{t}}_{{\mathrm{CO}}_2}^{-1}}.$$

(23)

By weighting the price per tonne of ${\mathrm{CO}}_2$ equivalent emission, contrails are considered in the multi-criteria trajectory optimization TOMATO. Following global approximations, as elaborated in the introduction, it is expected that the radiative impact of a single flight’s contrail exceeds those effects of ${\mathrm{CO}}_2$.

3. Results

3.1. Contrails in Trajectory Optimization

The most important effects of $R{F}_{\mathrm{Contrail}}$ are determined by the angle $\alpha$ between the length axis of the contrail and the incoming photons (see Figure 4), because it determines the travel distance of photons through the contrail. Furthermore, $R{F}_{\mathrm{Contrail}}$ depends on the order of magnitude of incoming radiation, depending on wavelength (solar > terrestrial), time of the day and year (summer, noon > winter, morning, and evening) and on latitude (equator > pole). Both effects are shown in Figure 9, because $\alpha$ and solar intensity change with daytime. Although the amount of backscattered power $P_{\downarrow b}$ dominates over the whole day, the sum of downward backscattered power $P_{\uparrow b}$ and absorbed power $P_{\downarrow \uparrow a}$ nearly compensates $P_{\downarrow b}$ between 7 a.m. and 5 p.m. Only during sunrise and sunset with horizontal photon transport (around 6 a.m. and 6 p.m.), is $R{F}_{\mathrm{Contrail}}$ significantly negative. The closer the contrail is to a pole (high latitudes), the longer the horizontal photon transport takes place. This is why, the ”sunrise–sunset effect” increases with increasing latitude, but the amount of incoming radiation (i.e., the amplitude in Figure 9 decreases with increasing latitude).

Although in Figure 9 the extinguished powers of only a single wavelength are shown, the statement of the results may be transferred to the narrow wavelength band around $0.2<\lambda <2\mathsf{\mu }$m. In Figure 10, simulations of the whole solar spectrum with significant contribution to the energy budget with wavelengths between $0.55\mathsf{\mu }\mathrm{m}\le \lambda \le 4.5\mathsf{\mu }\mathrm{m}$ (left) and of the most important terrestrial spectrum wavelengths between $4.5\mathsf{\mu }\mathrm{m}\le \lambda \le 21.5\mathsf{\mu }\mathrm{m}$ (right) are shown. Figure 10, left, indicates a decreasing solar effect with increasing wavelength due to a decreasing solar intensity and an increasing share of absorbed power. However, with increasing wavelength, positive values of $R{F}_{\mathrm{contrail}}$ become more dominant. Note, in Figure 9 and Figure 10, only the radiative properties of a single contrail with specific optical properties are calculated. The absorption efficiency $Q_a$, known as the possibility of an absorbing event within the contrail, strongly varies with ice-particle size and shape, which is why the fluctuation in Figure 10, left, is very contrail-specific.

The strong impact of ice-particle size on $Q_a$ causes an increasing terrestrial radiative forcing with increasing contrail lifetime (see Figure 2, Figure 3, and Figure 11, for increasing ice-particle radius with lifetime). Additionally, an increasing contrail width with lifetime causes a more distinct radiative forcing of older contrails. Note, the solar radiative forcing (Figure 11, left) is shown only for a single wavelength $\lambda =0.55\mathsf{\mu }\mathrm{m}$, whereas the terrestrial $R{F}_{\mathrm{Contrail}}$ is integrated over the whole significant terrestrial spectrum $4.5\mathsf{\mu }\mathrm{m}\le \lambda \le 21.5\mathsf{\mu }\mathrm{m}$. However, at $\lambda =0.55\mathsf{\mu }\mathrm{m}$ the solar intensities are maximum (see Figure 7). Hence, the main solar contribution to $R{F}_{\mathrm{Contrail}}$ is covered. Furthermore, the optical properties depending on shape, size, and number density of ice crystals and on wavelength are only parameterized for discrete wavelengths at $0.55\mathsf{\mu }\mathrm{m}$, $1.35\mathsf{\mu }\mathrm{m}$, $2.25\mathsf{\mu }\mathrm{m}$, $2.75\mathsf{\mu }\mathrm{m}$, $3.0125\mathsf{\mu }\mathrm{m}$, $3.775\mathsf{\mu }\mathrm{m}$, and $4.5\mathsf{\mu }\mathrm{m}$ [48]. These wavelengths are used as sampling points in the integration over the whole solar wavelength spectrum. these calculations have not to been carried out as functions of particle shape and size. That is why the solar radiative forcing (Figure 11, left) is only shown for a single wavelength. More detailed results are shown in [11].

Finally, the impact of the consideration of contrails in trajectory optimization is exemplified in Figure 12. Here, blue squares denote ice-supersaturated regions wherein contrail formation is very likely. For this example flight from Los Angeles (LAX) to Boston (BOS), different routes are calculated and compared with each other. The originally filed route (black, 4580 km ground distance) induced a contrail with a length of $L_{\mathrm{Contrail}}=1705$ km ground distance. The mean optical properties of this contrail have been estimated to $R{F}_{\mathrm{Contrail}}=0.305\mathrm{W}{\mathrm{m}}^{-1}$ which is $10.19\times {10}^{-9}\mathrm{W}{\mathrm{m}}^{-2}$ using Equation (17). Following Equation (23), this contrail must be weighted with $26.81\mathrm{t}{\mathrm{CO}}_2$ equivalent emissions. If this value is divided by the total length of the contrail ($L_{\mathrm{Contrail}}=1705$ km), we obtain an average amount of $0.0157\mathrm{t}{\mathrm{CO}}_2{\mathrm{km}}^{-1}$ equivalent emissions per km contrail. Considering $0.0157\mathrm{t}{\mathrm{CO}}_2{\mathrm{km}}^{-1}$ and €65 per tonne ${\mathrm{CO}}_2$ equivalent emission in a multi-criteria trajectory optimization, the contrail would have been completely avoided (red route in Figure 12 with 4730 km ground distance). A more harmonized trajectory with a reduced contrail length of $L_{\mathrm{Contrail}}=1374$ km ground distance constitutes the green route in Figure 12 with a total ground distance of 4246 km and a contrail weighting of $0.0125\mathrm{t}{\mathrm{CO}}_2{\mathrm{km}}^{-1}$.

3.2. Verification

To verify the contrail radiative forcing model, a comparison with the results of other publications is carried out whenever possible. Mostly, the optical thickness $\tau$ is used to describe the optical properties of a single contrail, measured vertically through the atmosphere. In the present study, the contrail is treated as a Gaussian plume without defined boundaries, but the optical thickness $\tau$ strongly depends on the contrail irradiated width. The larger the irradiated width, the smaller the optical thickness at constant irradiation, because the share of photons traveling through the contrail core is reduced. Furthermore, measurements of contrail optical properties are error-prone because of difficult measuring conditions at cruising altitude, where the contrail is hardly visible. For measuring the optical contrail properties, the shape, size, and number density of ice crystals have to be estimated. Therefore, the ice crystals must be captured and examined. Due to high aircraft speeds (≈$220\mathrm{m}{\mathrm{s}}^{-1}$), the crystals shatter when caught [58]. From this follows that only a few studies remain for a model verification. For example, Forster et al. [14] and Gounou et al. [13] published the corresponding contrail width to the optical thickness. This width can be interpreted as the irradiated width used in the present radiative forcing model. Hence, the weighted number ratios of extinguished to evaluated photons $S_{\mathrm{ext}}$, which is the sum of $S_{\downarrow a}$, $S_{\downarrow b}$, and $S_{\downarrow f}$ multiplied with the irradiated width, can be used for comparison with other publications. In

Table 1. Comparison of weighted number ratios $S_{\mathrm{ext}}$ for $\lambda =0.55\mathsf{\mu }\mathrm{m}$ of different ages of lifetime with Forster et al. [14]. The calculations are carried out for a mean turbulent environment with $\epsilon =5\times {10}^{-5}{\mathrm{m}}^2{\mathrm{s}}^{-3}$.

Study	Age [s]	Contrail Width/Horizontal Standard Deviation [m]	${\mathit{S}}_{\mathbf{ext}}$ [m]
Forster et al. [14]	1000	960	$223.6$
	2000	1680	$403.2$
	3500	2880	$720.0$
	6500	4320	$1036.8$
	$\mathrm{11,000}$	6000	$1020.0$
	$\mathrm{18,000}$	6480	$528.4$
present study	$\mathrm{18,645}$	$778.5$	$342.5$
	$\mathrm{35,798}$	$1078.5$	$453.0$
	$\mathrm{54,448}$	$1330.3$	$585.3$
	$\mathrm{72,011}$	$1539.9$	$677.5$
	$\mathrm{90,902}$	$1718.9$	$739.1$
	$\mathrm{120,224}$	$1976.7$	$889.5$

, values of the weighted number ratios $S_{\mathrm{ext}}$ of a non-sheared contrail in a mid-latitude summer atmosphere for several contrail ages derived by Forster et al. [14] are compared with the present weighted number ratios $S_{\mathrm{ext}}$ at $\lambda =0.55\mathsf{\mu }\mathrm{m}$. Although different atmosphere conditions are used by Forster et al., and despite the fact that Forster et al. integrated the radiative extinction over the whole solar and terrestrial spectrum, the results of the weighted number ratios are in the same order of magnitude.

The influence of the solar zenith angle on radiative extinction can be compared with Gounou et al. [13] as well as with Forster et al. [14]. Gounou et al. [13] published the radiative forcing $RF$ [$\mathrm{W}{\mathrm{m}}^{-2}$] of a dense blanket of contrails, integrated over the whole spectrum. Nevertheless, the shape of the extinction depending on solar zenith angle should be similar to the results observed in the present study. Figure 13 shows the typical minimum of radiative forcing at solar zenith angles around $\theta \approx 70$° as observed by Gounou et al. [13] and Forster et al. [14].

3.3. Discussion

In this study, a method is described to consider the environmental impact of individual condensation trails in a multi-criteria trajectory optimization. Thereby, three research questions have been answered:

• How can the conditions of contrail formation be estimated in a flight performance model?
• What is the radiative impact of the induced contrail per flight time step?
• How can the estimated contrail radiative impact be considered in a multi-criteria trajectory optimization?

Figure 1 can be used as a guideline for the implementation of the results of this study. Therein, the estimation of the conditions of contrail formation requires the atmospheric state parameters relative humidity and temperature (top left in Figure 1) as well as speed, thrust, and fuel flow from the flight performance module (top, center in Figure 1). A tendency towards aircraft speeds with higher fuel flows (resulting in a lower propulsion efficiency (Equation (2)), a lower slope of the mixing line G (Equation (1)) and therewith in a lower threshold temperature for contrail formation) has been identified to decrease the probability of contrail formation. High fuel flows may be induced by both too low and too high aircraft speeds. Since this decision (contrails formed or not) is crucial for further calculation steps, a highly resolved (in time and space) weather data set and a reliable trajectory calculation tool is indispensable.

For answering the second question, three modules are applied: The Contrail Life Cycle module (center left in Figure 1) for estimating the (micro-)physical properties, i.e., the dynamics of particle size distribution, particle shape distribution, contrail cross section, and sheared component during the entire lifetime, which again requires atmospheric data (top left in Figure 1) and the flight performance data (top, center in Figure 1). Beside aircraft type-specific and flight performance-dependent dimensions of the initial wake vortex, the atmospheric turbulence mainly influences the initial contrail dimensions (weak turbulence: small contrail) and the contrail lifetime (strong turbulence: short lifetime). Wind shear and vertical (up-)wind speed positively correlate with a longer contrail lifetime. The radiative properties are functions of the physical properties (center left in Figure 1) as well as of the atmospheric radiative transfer (top right in Figure 1), which quantifies the incoming radiation depending on solar zenith angle and Earth’s surface type. In general, the closer the sun is to the zenith, the higher the solar intensities irradiating the contrail. Furthermore, the warmer the Earth’s surface, the higher the terrestrial irradiances.

To consider multiple scattering events, a Monte Carlo simulation is set up in the Contrail Optical Properties module (center in Figure 1), where the position and the kind of extinction event are determined probabilistically.

Due to the separation of the contrail optical properties module (center in Figure 1) and the atmospheric radiative transfer (top right in Figure 1), the radiative extinction due to the contrail is calculated separately from the atmospheric radiative extinction. In the contrail radiative transfer module (lower right in Figure 1) the radiation, which has already been extinguished by the atmosphere itself without any contrail presence, is again extinguished according to the optical properties of the contrail. The advantage of this approach is that radiative properties of optically similar contrails (calculated in the “contrail optical properties” model) can be easily combined with different atmospheric conditions (estimated in the “atmospheric radiative transfer” model).

The third question is answered by transforming the contrail radiative impact into monetary values in order to compare the effect with other “costs” in the trajectory optimization. This has been done by transforming the radiative impact of the contrail into tonnes of CO${}_2$ equivalent emissions. The CO${}_2$ equivalent emissions are then transformed into monetary values by using the emission trading system. The weighting of the contrail impact in the trajectory optimization is adjusted with prices between €10 and €65 per tonne CO${}_2$ equivalent emission.

4. Conclusions

This study presents a method to granularly consider the ecological impact of individual condensation trails in trajectory optimization. If the implementation of the complete model as shown in Figure 1 exceeds the required calculation time, the following findings, which were drawn from the implementation into trajectory optimization, can be generalized: The identification of ice-supersaturated regions from weather data is very simple and can be implemented in the path-finding algorithm as layer with additional costs to simply avoid contrail formation. If the costs for the detour around ice-supersaturated regions exceed the environmental costs of the contrail, the following costs can easily be included in the path finding. The contrail lifetime and therewith its radiative impact is reduced in regions with slightly higher turbulence (with eddy dissipation rates $\epsilon \approx {10}^{-4}{\mathrm{m}}^2{\mathrm{s}}^{-3}$, which still corresponds to light turbulences according to ICAO categorization [59]). Furthermore, avoiding regions with strong upwinds, which typically occur in cumulus clouds or thunderclouds, reduces the contrail lifetime, because strong upwinds lift the contrail in cold and ice-supersaturated regions.

In mid-latitudes, even during the daytime, the warming effect of the contrail dominates for two reasons: A significant amount of absorbed terrestrial powers and a strong forward scattering of solar powers reduces the amount of backscattered solar powers. During sunrise and sunset, when solar radiation is transferred horizontally through the contrail, the cooling effect often dominates, especially if the contrail longitudinal axis is parallel to the direction of the solar photons. This significant impact of the track angle (i.e., the contrail length axis) on the radiative properties determines the travel distance of horizontal photons through the contrail. Although suggestions of preferring north–south-oriented routes to east–west routes are less helpful in daily operations, the effect is considered in the costs of the pathfinding algorithm. To minimize the warming effect, flying above cold surfaces, e.g., snow, reduces the amount of terrestrial powers which might be absorbed by the contrail. This effect is also considered in the trajectory optimization.

Therewith, a guideline for a precise quantification of the contrail radiative impact is provided. Additionally, generalized findings of this approach for a simple consideration of the contrail impact in an arbitrary pathfinding algorithm for aircraft trajectory optimization has been derived. Both will be used in future applications to improve the quantification of environmental costs in aircraft trajectory optimization tools.