Aircraft routing with minimal climate impact : the REACT 4 C climate cost function modelling approach ( V 1 . 0 )

In addition to CO2, the climate impact of aviation is strongly influenced by non-CO 2 emissions, such as nitrogen oxides, influencing ozone and methane, and water vapour, which can lead to the formation of persistent contrails in ice-supersaturated regions. Because these non-CO 2 emission effects are characterised by a short lifetime, their climate impact largely depends on emission location and time; that is to say, emissions in certain locations (or times) can lead to a greater climate impact (even on the global average) than the same emission in other locations (or times). Avoiding these climate-sensitive regions might thus be beneficial to climate. Here, we describe a modelling chain for investigating this climate impact mitigation option. This modelling chain forms a multi-step modelling approach, starting with the simulation of the fate of emissions released at a certain location and time (time-region grid points). This is performed with the chemistry–climate model EMAC, extended via the two submodels AIRTRAC (V1.0) and CONTRAIL (V1.0), which describe the contribution of emissions to the composition of the atmosphere and to contrail formation, respectively. The impact of emissions from the large number of time-region grid points is efficiently calculated by applying a Lagrangian scheme. EMAC also includes the calculation of radiative impacts, which are, in a second step, the input to climate metric formulas describing the global climate impact of the emission at each time-region grid point. The result of the modelling chain comprises a four-dimensional data set in space and time, which we call c imate cost functions and which describes the global climate impact of an emission at each grid point and each point in time. In a third step, these climate cost functions are used in an air traffic simulator (SAAM) coupled to an emission tool (AEM) to optimise aircraft trajectories for the North Atlantic region. Here, we describe the details of this new modelling approach and show some example results. A number of sensitivity analyses are performed to motivate the settings of individual parameters. A stepwise sanity check of the results of the modelling chain is undertaken to demonstrate the plausibility of the climate cost functions.


Introduction
The anthropogenic origin of a substantial contribution to observed climate change is well established (e.g.IPCC, 2007).The challenge is how to deal with climate change and to find and evaluate mitigation strategies.Air traffic has a significant contribution to total anthropogenic climate change (IPCC, 1999;Berntsen and Fuglestvedt, 2008;Lee et al., 2010;Burkhardt and Kärcher, 2011) and a significant part of its contribution arises from non-CO 2 emissions, e.g. from changes in ozone, methane, cloudiness and others.
Published by Copernicus Publications on behalf of the European Geosciences Union.
These non-CO 2 effects are characterised by a high temporal and spatial variability; that is, their impact on climate depends not only on the amount of emitted species, as in the case of CO 2 , but also on the time and region where the emissions take place.The formation and persistence of contrails depend on both aircraft and fuel parameters and meteorological conditions (Schumann, 1996), such as ice supersaturation.These ice-supersaturated regions are locally and temporarily very confined and show a large variability (Spichtinger et al., 2003;Gierens and Spichtinger, 2000) and a dependence on the prevailing weather conditions (Irvine et al., 2012).Contrails may persist for a long time under favourable weather conditions (Minnis et al., 1998;Burkhardt and Kärcher, 2009).The variability of the climate impact from NO x emissions with respect to weather conditions has not been investigated yet.Climatological studies (Grewe and Stenke, 2008;Köhler et al., 2008;Frömming et al., 2012) show a distinct altitude and latitude variability.For an individual weather situation this variability is probably enlarged for all non-CO 2 effects.For example, NO x emissions which take place in a region with cloud formation and rain will have a significantly lower impact than an emission taking place in a region with upwelling, increasing the lifetime of the emitted species.
Within the EU project REACT4C (http://www.react4c.eu/;for abbreviations see Table A1), we quantify this variability and use it to develop possible strategies whereby aircraft are routed to minimise their total climate impact.It is likely that the fuel consumption increases for these aircraft trajectories, because in general (although not always) aircraft currently take routes that are close to the minimum fuel (and hence minimum CO 2 emission) route.In the case of climateoptimal routing, the gain from non-CO 2 effects counteracts the CO 2 induced warming, at least to some extent.How strong the compensating gain from the avoidance of non-CO 2 climate effects needs to be also depends on the objective.It will need to be stronger if the objective is on long-term climate change, since CO 2 effects are then more pronounced.Thereby the choice of climate-optimal route depends on the chosen time frame (or time horizon) and the adopted indicator of climate change.
Basically this optimisation is based on two major steps, i.e. the calculation of 1. climate cost functions and 2. aircraft trajectories optimised on the basis of the climate cost functions, where the climate cost functions (CCF) are specific climate metrics, i.e. climate impacts per unit emission.The idea of weather-specific re-routing of air traffic for the benefit of climate has been addressed before (Sausen et al., 1994;Mannstein et al., 2005;Schumann et al., 2011;Sridhar et al., 2012).However, none of these studies included such a broad range of effects, as addressed in this study: contrails, carbon dioxide, ozone, methane, ozone from methane changes, and water vapour.The changes in ozone arising from changes of its precursor methane are also called primary-mode ozone (PMO).
The idea of REACT4C (Matthes, 2012;Matthes et al., 2012) is first to concentrate on frequently occurring daily weather patterns (Irvine et al., 2013), for which a detailed analysis is performed with a new modelling approach, which we describe here in detail.Results will be published in companion papers.The methodology is outlined in Figs. 1 and 2. As a first step, we are concentrating on the North Atlantic region, including most of Europe and North America.For this region a weather classification is performed taking into account specific air traffic routes (Irvine et al., 2013).For each type of weather pattern a representative weather pattern is selected.For this 1-day weather pattern, a time-region grid (Fig. 3) is defined and for each time-region grid point the global climate impact of an emission is calculated for different emission times.The time-region grid covers the flight tracks over the North Atlantic and the main cruise levels.
The climate cost functions are calculated with the chemistry-climate model EMAC (version 2.42), which additionally includes two important submodels: ATTILA, a Lagrangian transport scheme (Reithmeier and Sausen, 2002) and AIRTRAC (version 1.0, Frömming et al., 2013; see Supplement), which calculates contributions from additional emissions to concentrations based on ATTILA.The EMAC submodel AIRTRAC is specifically developed for this study and described in Sect.3. The determination of the climate cost functions (Fig. 1, highlighted box in the left column) includes first the calculation of the contributions of additional emissions to atmospheric concentrations (nitrogen oxides, ozone, methane, contrails, water vapour, carbon dioxide) and contrail properties, and second the calculation of the radiative impact over a time period of weeks leading to an approximate annual mean instantaneous radiative forcing.Third, we use a correlation between instantaneous and adjusted radiative forcing to obtain the latter as a more reliable basis for the expected climate change (Sect. 3.4).This is eventually used as input to climate response formulas to obtain a set of metrics per unit emission, i.e. climate cost functions, which is interpolated back to the original EMAC grid.
In Sect.3.5 we relate individual metrics to political objectives and optimisation problems.The climate cost functions are used in the next step (Fig. 1) by a flight planning tool (SAAM) to obtain aircraft trajectories and respective emissions as well as the reduction in the climate impact due to the operational changes in aircraft trajectories.The climate cost function approach aims at reducing the contribution of air traffic on climate.The results for individual weather patterns can be multiplied by their frequencies to obtain an estimate of the total climate impact reduction as a result of the REACT4C re-routing strategy.Since this leads to changes in the background concentration and production efficiencies, e.g.ozone production per NO x molecule, other emitters such (not described in detail) Fig. 1.Illustration of the model chain.Models are given in the ovals.Data and definitions are given in rectangles.Climate cost functions are highlighted (dark blue), as are major results (green).The light-blue box at the bottom indicates a possible extension of the modelling chain, which is not covered here.
as industry or traffic might have a larger climate impact (Grewe et al., 2012a).This can be investigated by applying atmospheric-chemistry models to obtain the total mitigation gain (Fig. 1, green boxes).This is planned in the project RE-ACT4C, but is not part of this publication (light-blue box in Fig. 1).The whole methodology is based on operational models which were extended.The operational models are briefly introduced in Sect.2, and the new modelling approaches are described in detail in Sect. 3 (climate cost functions).In Sect. 4 we provide a comparison to other studies and a sanity check of our modelling approach.

Base models
The REACT4C modelling approach is based on a number of models which have been applied many times previously.These models are combined with new approaches.Therefore we only briefly describe the base models and concentrate on a detailed description of the new approaches in Sect.3.

Atmosphere: EMAC
The ECHAM/MESSy Atmospheric Chemistry (EMAC) model (here version 2.42) is a numerical chemistry and climate simulation system that includes submodels describing tropospheric to middle-atmospheric processes and their  interaction with oceans, land and human influences (Jöckel et al., 2006).Here we use the second version of the Modular Earth submodel System (MESSy2) to link the individual physical and chemical processes described in submodels (Jöckel et al., 2010).The MESSy interface also provides standard and well-tested routines, such as those for data extraction and grid remapping, which facilitates the implementation of new submodels.We run the model in a T42L41 resolution, which is approximately 2.8 • × 2.8 • latitude-longitude resolution and 41 layers in the vertical from the surface to an upper layer centred at 5 hPa.
The core model for the atmosphere is the fifth generation European Centre Hamburg general circulation model ECHAM5 (Roeckner et al., 2006).The chemistry is described by the submodel MECCA (version 3.2) (Sander et al., 2011).More detailed information, including references, about the model system is available from http://www.messy-interface.org.

Aircraft routing and emissions
The simulation of the flow of air traffic is performed with the System for traffic Assignment and Analysis at a Macroscopic level (SAAM) to which the Advanced Emission Model (AEM) is coupled.An overview on the simulation system SAAM and AEM is given in Fig. 2 and discussed in more detail in Sect.2.2.3.
It allows for the creation, change and design of air traffic route networks with their possible associated constraints (e.g.restrictions from the route availability document -RAD or flight level constraints).From any traffic demand (basically worldwide airport origins and destinations, aircraft types and departure times), a set of full 4-D aircraft trajectories (space and time) is generated.
The best choice of these 4-D aircraft trajectories is made using optimisation with an objective function minimising a mathematical cost that can be based on either economical values derived from flight time and aircraft operating cost or on the climate impact (see below).Constraints concerning conflict avoidance between all 4-D aircraft trajectories can be switched on or off in the optimisation model.Climate impact minimisation has not been performed before with SAAM, but other optimisation problems were considered, for instance a balancing of the air traffic controller's work load (Champougny et al., 2001).
Aircraft performance included in SAAM uses Base of Aircraft Data (BADA) version 3.9 (see http://www.eurocontrol.int/services/bada), which also provides fuel flow based on aircraft-engine characteristics and assumptions like mean aircraft load and weight.

AEM
The advanced emission model (AEM 2.5.0) has been developed to estimate the mass of fuel burnt, and emissions produced by a specific aircraft-engine configuration for a specific 4-D aircraft trajectory (Eurocontrol, 2013).Emissions are calculated for CO 2 , H 2 O, NO x , SO x , CO and unburned hydrocarbons (HC), e.g.benzene (see also http: //www.eurocontrol.int/services/advanced-emission-model).

Application of SAAM with AEM
Figure 2 shows, in more detail, the information flow and coupling of SAAM and AEM.First, flights are selected from a database -including the city pair connection, departure time and aircraft and a reference aircraft trajectory calculated -based on minimum costs.Second, routes are generated randomly.Here we apply a procedure which generates alternative routes by blocking of air space areas in addition to a randomly changed cruise altitude.Three different grid sizes (Fig. 4) for the blocking of the air space were tested.The SAAM standard grid size (grid 0) leads to significant deviations from the great circle, whereas the smaller cells only provide small deviations.Hence, the resolution with the larger grid cells (grid 0) is better suited for the generation of alternative routes and taken for the optimisation.This procedure leads to a randomly chosen variation of 16 additional routes in the horizontal with 4 additional options for the cruise altitude (3 below the original cruise level and 1 above) and hence 84 (17×5−1) alternative routings for each reference flight (i.e.city pair connection for a specific time).The number of the blocked grid cells has been varied to test the sensitivity of this parameter to the optimisation process.We found that an increase in the number of blocked cells to 18 only shows a minor change in the optimal solution with differences well below 1 %.
The third step (Fig. 2) is the calculation of the 4-D aircraft trajectory, which includes a performance calculation and takes into account wind fields.This is followed by the fourth step, namely the calculation of the emissions by using the AEM model.This leads to a large set of aircraft trajectories.Based on these data, an optimal air traffic flow is determined.The optimisation is either done with respect to economic or climate costs (climate impact).In the case of ignoring any conflicts, the optimisation simply chooses the minimum among the alternative routes for every reference flight (i.e.city pair connection at a specific time).In the case of conflict avoidance, the routes depend on each other and a linear programming is applied: with the flight index i (and k, respectively): 1 ≤ i ≤ n (n is the number of city pair connections), j (l, respectively) the route option index (≤ 85; see above), the precalculated conflict matrix CM ≥ 0 (i.e.number of conflicts) and the costs C i,j .These costs are either defined for economic optimisation or climate optimisation.For economic optimisation the costs are with F i,j and T i,j the fuel consumption and flight time on route option (i, j ) in [kg] and [min], respectively, and C fuel = EUR 0.75 kg −1 , C time = EUR 25 min −1 .For climate optimisation the costs are defined as where the aircraft trajectory of route option (i, j ) is divided into flight legs with index m, which is of the order of 10

Methodology: climate cost functions
In this section we describe the calculation of the climate cost functions, starting with the Lagrangian approach and then followed by the definition of the time-region grid points and the chemistry and microphysics of the atmospheric processes, radiation changes and the climate impact calculations.A summary is presented in Sect.3.6 and Table 1.The atmospheric process modelling is performed in the EMAC submodels AIRTRAC and CONTRAIL.Further documentation of the program structure, subroutines and namelists is given in the Supplement.

ATTILA properties
time-region 2 time-region 3 time-region 4 time-region 5 time-region 1 Fig. 5. Properties of an individual Lagrangian air parcel (trajectory) in the submodel AIRTRAC.The possibility to calculate numerous time-region grid point emissions in one simulation is a key aspect for computationally efficient coding.

A Lagrangian approach: trade-off between computational efficiency and process accuracy
The calculation of climate cost functions requires a calculation for each of the predefined time-region grid points.Hence, the desired resolution of the climate cost function in time and space determines the number of climate cost function calculations.We denote the number of time-region grid points by N TR and denote the individual time-region grid points by TR i (i = 1, . . ., N TR ).We chose a Lagrangian approach since it allows the inclusion of a multitude of cost function calculations in a single EMAC simulation.Each air parcel trajectory is characterised by its position at any time and includes an arbitrary number of properties P .We assign for each time-region grid points TR i a set of n properties, i.e.P ((i − 1) • n + 1), . . ., P (i • n). Figure 5 shows examples for five time-region grid points (each coloured differently) a set of eight properties, namely the contribution of the emissions to the chemical species NO, NO 2 , HNO 3 , O 3 , CH 4 , and H 2 O and contrail coverage (see Supplement for a complete list of the 13 (nlgtrac= 13) properties).
The Lagrangian approach is based on detailed modelling of the background processes within EMAC and an additional, and to some degree simplified, simulation of the contributions from emissions taking place in the respective time regions.Information from the detailed modelling is transferred to the air parcel trajectories and the contributions are calculated based on the results from the detailed process modelling within EMAC.
In principle two modelling approaches are applicable: (1) the Lagrangian approach introduced above and (2) a perturbation approach, which includes a base case simulation as well as an additional full simulation with the base model (here EMAC) for each time region, including additional emissions from the time region.Here we chose the Lagrangian approach since it better meets our objective: obviously the Lagrangian approach has the advantage of being numerically efficient, since many time regions are calculated  3 and 4) Saturation effects for NO x chemistry, contrail cirrus and water vapour are less important than the 3-D structure of the response to the emission This verification will be performed in a followup publication, since it includes an analysis of the implication on air traffic emissions, which is beyond the scope of this work (see also shaded box in Fig. 1).

Chemistry
Tagging Calculation of contribution of the time-region emission to the atmospheric composition independently of each other (no feedback to background chemistry) Linearisation of the specific reaction rates for small perturbations (Eq.9) is valid The overall objective of our modelling approach is to minimise the contribution of air traffic to climate change.It is important to stress that the contribution calculation (Grewe et al., 2010) is much better suited to address this objective than the perturbation method, since it does not lead to misinterpretations of the results due to non-linear compensation effects in the chemistry (Grewe et al., 2012a).
Changes from the time-region emissions do not feed back to the base model processes.This ensures an identical background meteorology and chemistry for all time-region simulations; that is, we perform quasi-chemical-transport model (QCTM) simulations (Deckert et al., 2011).
The second approach, perturbation simulations for each time region, would have the advantage of having the same degree of detail in the description of background processes and the perturbation.However, it is extremely time-and resourceconsuming, and only its total impact is calculated, not the contribution of the time region (see above and also Grewe et al., 2010Grewe et al., , 2012a;;Grewe, 2013).

Definition
Figure 3 shows the various grids used in our approach.The EMAC grid has a resolution of approximately 2.8 • longitude × 2.8 • latitude (black).On this grid we overlay the timeregion grid (red triangles), covering the North Atlantic region for the area where the optimisation with SAAM is performed (see Sect. 2.2.3).It consists of six longitudes, seven latitudes and four pressure levels for three points in time (see Table 2 for details).This comprises 504 grid points, meaning 504 calculations for the impact of emissions on the atmosphere have to be performed.We take into account water vapour emissions and nitrogen oxide emissions (see Table 3).The emissions are released for one time step (15 min) in the EMAC gridbox, in which the respective time-region grid point is located.Additionally, we consider properties of conventional present-day aircraft and fuel, such as the aircraft's overall propulsion efficiency, the combustion heat and the water vapour emission index (see Table 3).
For each time-region grid point, emissions are partitioned on 50 air parcel trajectories (see also inlay in Fig. 3), which are randomly distributed in the EMAC grid box in which the time-region grid point is located.Technically this is handled via the submodel TREXP and controlled by a namelist (Jöckel et al., 2010).Then the impact of emissions is calculated by the AIRTRAC submodel (box model on the air parcel trajectories), the radiative forcing calculated and the climate impact (measured by climate metrics; see also Fig. 1) is associated with the time-region grid point.Finally these climate impact data on the time-region grid are interpolated to the EMAC grid, which is then the final climate cost function grid (Fig. 3).

Sensitivity to the number of air parcel trajectories
Here we investigate the number of air parcel trajectories at each emission point that are required in order to obtain unambiguous results.If the number of air parcel trajectories is too low, then only a few possible paths for the air parcels carrying the emission will be considered, and thus the results can become very noisy and random.
For this sensitivity study, 50 air parcel trajectories were released for 24 of the time-region grid points, located at 200 hPa at 45 • N and 50 • N and at 400 hPa at 40 • N and 35 • N.For each of these time-region grid points, the mean NO x (NO + NO 2 ) mixing ratio over the first month of integration (i.e.January) was calculated over all 50 air parcel trajectories.We tested how strongly the results differ when less (between 2 and 48) air parcel trajectories were released.For each number of air parcel trajectories (between 2 and 48) 100 sub-samples out of the total 50 air parcel trajectories were randomly generated.For each sub-sample, the mean NO x mixing ratio was calculated and compared to the mean over all 50 air parcel trajectories.The relative deviations for these sub-samples as a function of number of air parcel trajectories is shown for one time-region grid point in Fig. 6 (top, left).Mean NO x values differ by up to 40 % from the result for all 50 air parcel trajectories if only 2 air parcel trajectories are used.With an increasing number of air parcel trajectories, the deviations decrease and converge towards the mean value over 50 air parcel trajectories.The same diagnostic is shown for the corresponding mean ozone mixing ratios (Fig. 6, upper right).This diagnostic shows a similar behaviour to NO x , but with a smaller initial spread in the results for small numbers of air parcel trajectories.In the bottom panel of Fig. 6, the standard deviation as a function of the number of air parcel trajectories shows a sharp decrease of the NO x and ozone standard deviation from around 20 % to 5 % and 12 % to 4 %, respectively, up to an air parcel trajectory number of around 20.This figure confirms that the deviation in the results of ozone is smaller than in the NO x values also for the mean over the 24 time-region grid points at different locations.
While the mean standard deviation of the results is almost 20 % for NO x and 12 % for ozone when using only 2 air parcel trajectories, it drops quickly with an increasing number of air parcel trajectories and both the mean values and the extremes lie below 10 % for more than 20 air parcel trajectories.Thus, this analysis suggests that the potential error in the results for NO x and ozone almost converged for 50 air parcel trajectories and is lower than 10 % if 20 or more air parcel trajectories are used.

Sensitivity to temporal and horizontal resolution
We have tested the impact of the temporal and horizontal resolution of the time-region grid on the climate cost functions.
For the test of the temporal resolution, we included an emission time at 09:00 UTC in addition to the standard emission times at 06:00 UTC and 12:00 UTC at 250 hPa and between 45 • N and 50 • N. The interpolated values are compared to the simulated 09:00 UTC values.The NO x and ozone masses vary by roughly ±40 % and ±25 %.
Concerning the horizontal resolution we have added emission locations in between the chosen grid, i.e. a latitudinal shift of 2.5 • and longitudinal shift of 7.5 • at 250 hPa and 12:00 UTC.The results show a variability of the order of 50 % and 35 % for NO x and ozone, respectively.
Hence the horizontal interpolation is more critical than the temporal interpolation.The intercomparison further shows that the resolution of the climate cost function is crucial to the optimisation of the air traffic system with respect to its climate impact.However, the variation of the climate impact from NO x emissions varies by one order of magnitude and is hence larger than the possible interpolation error.

Chemistry: NO x , O 3 (+PMO), OH, CH 4
We consider a NO x emission of 5 × 10 5 kg for each timeregion grid point.This influences ozone production and loss as well as the HO x partitioning and hence methane loss.
From the detailed modelling in EMAC, we obtain relevant production (P) and loss (L) terms.Those terms are then used in the following tagging approach.A prerequisite of this tagging approach is a complete and disjunct partitioning of the emissions into categories.Here we consider two emission categories, a background (b) and the additional emission (e) for calculation of the climate cost function.For this emission category e the contribution of this category to the ozone production via the reaction is calculated (according to Grewe et al., 2010): where the superscripts b and e indicate background values and values specific for the tagged emission category, respectively.P O 3 denotes the ozone production rate in [mol mol −1 s −1 ].All species are given in mixing ratios [mol mol −1 ].The reaction rate for loss terms of tagged species is determined in the same manner, e.g. for the reaction the loss term L e O 3 is Note that the first term in brackets includes the loss of background ozone by NO e 2 and that Eqs. ( 7) and ( 8) only represent the contributions from an emission e for the Reactions (R1) and (R2) without any further assumptions, e.g.without any linearisation of chemical processes.
For our approach we now introduce simplifications.First we combine background nitrogen species (all N-compounds except NO, NO 2 and HNO 3 ) to a family and calculate the contributions of time-region emissions according to this family concept, with the assumption that the emissions are small enough that the specific reaction rates are unchanged, i.e. driven by two different chemical families: one based on NO x , the other taking into account all other loss processes, which leads to the following differential equation for ozone: This implicitly includes the assumption that the NO-to-NO 2 ratio of the NO x emitted in the time region is equal to the ratio of the background NO x .After NO x is emitted, we take into account an exchange with HNO e 3 , which will eventually be washed out.
Figure 7 shows the evolution of the contributions to the atmospheric burden [kg] of various species as an example for an emission at 200 hPa, 30 • N and 75 • W. Nitrogen oxides are completely washed out within a month.Ozone increases as long as enough NO x molecules are available and is subsequently destroyed when this is no longer the case.An analogous approach is used for methane.We take into account the most relevant reactions regarding the concentration of OH and HO 2 : Production of OH: Loss of OH: Production of HO 2 in addition to Reactions (R6) and (R7): Loss of HO 2 in addition to Reactions (R4), ( R5) and (R10): We are not considering aircraft contributions for H 2 O, CO, RH and CH 4 for this approach, since their effects on OH are considered to be small.The production and loss of OH e , i.e.OH of the respective emission category, follow in analogy to Eq. ( 7): with x ( 24) x ( 27) The methane loss caused by the contribution of aircraft emissions to the OH concentration is then Figure 7 shows, in addition to NO x and O 3 , the evolution for methane.Methane decreases first because of an increase in OH via Reaction (R1).When NO x is removed from the atmosphere, the increase in ozone concentrations leads to an enhanced methane reduction.

Aircraft-induced cloudiness
Contrails form in the atmosphere when the ambient air at flight levels is sufficiently cold and moist (Schmidt-Appleman criterion, SAC ;Schumann, 1996).Once formed, contrails may persist if the air is supersaturated relative to ice and evolve into contrail cirrus, i.e. they lose their line-shaped structure.Here we generally refer to contrails and do not distinguish between line-shaped contrails and contrail cirrus.We determine the atmospheric ability to form persistent contrails, i.e. the potential contrail coverage, instantaneously within the climate model at each time step following the approach of Burkhardt et al. (2008) and Burkhardt and Kärcher (2009).The potential contrail coverage, b co , is the fraction of an EMAC grid box which can be maximally covered by contrails.It is calculated as the difference between the maximum possible coverage of both contrails and cirrus (b co+ci ), and the coverage of natural cirrus alone (b ci ): Here, r denotes the EMAC grid mean relative humidity.r ci and r co are critical relative humidities above which a fraction of the EMAC grid box is covered by cirrus and is ice-supersaturated, respectively; r sat = 1 is the relative humidity at saturation.The relative humidity r * = r sat − (r ci − r co ) 2 /(r sat − r ci ) (see Burkhardt et al., 2008, and supplementary material therein).
The critical relative humidity for contrail formation, r co , is calculated via with r SAC being the relative humidity over ice at which the Schmidt-Appleman criterion is fulfilled during the mixing of aircraft exhaust gases and ambient air, and r nuc being the homogeneous freezing threshold.As contrails often form prior to the formation of natural cirrus a = 0.9 is chosen (following Burkhardt et al., 2008).
The potential contrail coverage is transferred to the air parcel trajectories.Then the actual contrail coverage is calculated depending on whether air traffic is actually taking place in the respective EMAC grid box.Further physical processes such as contrail spreading, sedimentation of ice particles, water uptake and sublimation are parameterised.Contrails are described by their coverage (b) and water mixing ratio (m).
The prognostic equation for the contrail coverage is the sum of newly formed contrails and the spreading of preexisting ones: with the following boundary conditions that the resulting contrail coverage b fulfils: The newly formed contrails cover an area according to the initial contrail dimensions, such as the width W 0 and the length of the contrail in the EMAC grid box (D; see also Table 4).We assume a flown distance A, and hence an initial contrail length of ( 38) where A is the gridbox area.The resulting new contrail coverage tendency is then The spreading of the contrails is parameterised according to Burkhardt and Kärcher (2009) depending on the vertical wind shear: The prognostic equation for the contrail ice mass mixing ratio m includes the formation of new contrails, sedimentation (or precipitation) of ice mass, deposition of water vapour on the contrail ice particles, and sublimation (Burkhardt and Kärcher, 2009): Similar to Ponater et al. (2002) we assume that the newly formed contrail ice water mixing ratio depends on the condensation rate in the contrail-covered part of the EMAC grid box c co , which is defined analogously to the condensation rate for natural clouds c cl : Extending over Ponater et al. (2002), however, contrail ice is also subject to physical sinks: the sedimentation rate of ice particles in the contrail is parameterised according to Heymsfield and Donner (1990) based on the vertical divergence of the flux of ice particles F in [kg m −2 s −1 ]: where ν is the falling velocity [m s −1 ] of ice particles, which are parameterised with α = 3.29 and β = 0.16, and ρ is the air density in [kg m −3 ].
The sublimation and growth is parameterised according to a relative decrease or an increase of the potential contrail coverage, resulting from a decrease or increase of the relative humidity, respectively: Note that the potential contrail coverage is available at a higher resolution than the RF calculated on the time-region grid points.This information is used after the mapping from the lower resolution of the time-region grid to the higher resolution EMAC and CCF grid (Fig. 3) to ensure that contrail climate cost functions differ from zero only where contrails form in EMAC.

H 2 O and CO 2
For every time region, a water vapour emission is taken into account (Table 3).This emitted water vapour (H 2 O in [mol mol −1 ]) is transported via Lagrangian transport, like the other tracers.Only loss processes are considered: where pr is the precipitation rate, i.e. the water vapour loss in [mol mol −1 s −1 ] in the respective EMAC gridbox due to rain and snowfall, and H 2 O tot is the respective EMAC grid box total water vapour in [mol mol −1 ].
Carbon dioxide emissions are assumed to lead to a wellmixed enhancement of the carbon dioxide mixing ratio because of its long perturbation life time.The temporal evolution of the concentration change (decay) is given for a unit of fuel used CO 2 (t) in [kg(CO 2 ) (kg(fuel)) −1 ] following Fuglestvedt et al. (2010) and Forster et al. (2007): (see Table 5 for details).

Radiative forcing
Table 6 gives an overview on the RF calculation for the individual species.The RF calculation for CO 2 is based on Fuglestvedt et al. (2010) and includes a simple linearised conversion factor between the change in its atmospheric mass and the RF of 1.82×10 −12 mW m −2 (kg(CO 2 )) −1 .The contrail RF is derived from the global mean radiation flux changes F at the tropopause: where T =1 yr.
Ozone RF is calculated analogously, except that the instantaneous RF is converted into the adjusted RF (details are described in the Appendix) via an analytical formula: The unitless functions f 1 and f 2 describe the relation between RF adj and RF inst for various times of the year (f 1 , seasonal cycle) and various perturbation altitudes (f 2 ) and are derived from additional idealised simulations.
Methane RF is derived from the methane mass changes calculated explicitly (Sect.3.3.1)by applying the respective IPCC formula (Shine et al., 1990).PMO radiative forcing is derived from the methane radiative forcing by a constant factor of 0.29 (Dahlmann, 2012).Water-vapour-adjusted RF is calculated based on results from Grewe and Stenke (2008, see Fig. 8).Grewe and Stenke (2008) investigated the consequences of sustained water vapour emissions at different atmospheric locations between the surface and 50 hPa and between the North Pole and the tropics.The results show that the adjusted RF depends on the lifetime of the perturbation and a linear relationship between the resulting atmospheric mass change and RF (Fig. 8).Note that Fig. 8 does not show a direct correlation between the water vapour mass changes and the RF, but does show that a given emission leads to a change in both the atmospheric mass and a RF, which both depend on the location where the emission takes place.Here, we have the same assumption, which is a constant emission at different locations (i.e.time-region grid points) and calculate the mass change explicitly (Sect.3.3.3)so that we can use the mass-RF relationship (Fig. 8).

Climate metrics and emission scenarios for the REACT4C objective
We are aiming at minimising the air traffic's climate impact by alternative routings.However, the wording "climate impact" or "climate change" is not well defined and incorporates a variety of possible interpretations.Hence, we clarify our objectives and from that derive adequate climate metrics and emission scenarios (Grewe and Dahlmann, 2012).
There is no uniquely applicable climate metric, so we implement several possible choices (climate objectives) in order to examine the implications on routing for each of them.
Our approach differs only slightly from previous considerations, which stress the importance of choices on the emissions, metrics and time horizons (see e.g.Shine et al., 2005;Fuglestvedt et al., 2010).Here, we move the focus slightly: we spell out three political questions and group appropriate combinations of emissions, metrics and time horizons.Still, choices have to be made.The emissions, metrics and time horizons are those which are frequently used: pulse (P), sustained (S) and future emission scenarios (F); absolute global warming potential (AGWP); absolute global temperature potential (AGTP); and average temperature response (ATR) for time horizons of 20, 50, and 100 yr.We aim at reducing the climate impact of air traffic by introducing a new routing strategy (more details in Sect.A2).That implies that this strategy is applied everyday for air traffic.We interpret climate change as the global mean temperature change.The definition of the strategy and the focus on global mean temperature change as a climate change indicator implies an emission scenario and a climate metric (see Table 7): the first question (Q1) is "what is the short-term climate impact of the REACT4C re-routing strategy?".The appropriate emission scenario is a best estimate for the future air traffic, and the average temperature response (ATR) for a 20 yr time horizon, H = 20, is a suitable climate indicator, which is abbreviated as F-ATR20 (Future emission-scenariobased Averaged Temperature Response on a 20 yr time horizon): Since estimates of the future air traffic are naturally uncertain, one can argue that a sustained emission is an adequate assumption (abbreviated as S-ATR20).A third option is to replace the mean temperature change over the time horizon H with the temperature change at the time horizon H , i.e. the absolute global temperature potential for sustained emissions (S-AGTP20) (see Supplement and Fuglestvedt et al., 2010).Furthermore, as a last option, with the assumption that the climate sensitivities are equal to 1 K (W m −2 ) −1 for all species (see Eq. A1), the absolute global temperature potential for the time horizon H = 20 can be approximated by the integrated radiative forcing for a pulse emission, i.e.P-AGWP20.
The second objective, i.e. the second question (Q2) is "what is the long-term climate impact of the REACT4C rerouting strategy?".This changes the focus from the near future to longer term effects, and hence, except for the time horizon (H = 100), the emission scenarios and metrics are the same as for Q1.
For reasons of completeness, we add here the third question (Q3): "what is the medium-range climate impact of a present-day REACT4C re-routing decision?".Here, a medium (H = 50) time horizon of 50 yr is addressed, but The choice of the objective has a consequence on the importance of individual atmospheric processes (Table 7).From question Q1 to Q3 the focus shifts from short-to long-term effects and hence from contrails and ozone impacts to CO 2 only.In order to restrict future analysis to both a minimum of metrics and a maximum in the range of uncertainties from the metrics, we suggest focusing on the metrics F-ATR20, F-ATR100, P-GWP20 and P-GWP100 in future applications.

Summary on climate cost function calculation
We have set up a modelling approach linking potential emissions at locally and temporarily confined regions to their climate impact, measured with climate metrics.This procedure follows previous approaches (e.g.IPCC, 1999;Grewe et al., 2007;Grewe and Stenke, 2008;Lee et al., 2010); however, it differs in some details, taking into account new findings.Table 1 shows the main features of this modelling chain.We used different grids (see also Fig. 3), starting from the grid of the base model EMAC, in which we included 504 selected time-region grid points which cover the North-Atlantic region.From the area immediately around each of the timeregion grid points, we started 50 air parcel trajectories each and allocated the resulting climate impact to this grid point.An interpolation to the original EMAC grid provides the final climate cost function grid.The impact of the interpolations and definition of the time-region grid is tested in Sect.3.2 (see also Table 1).
Physical and chemical processes are calculated on the air parcel trajectories by extracting process information from the detailed EMAC model.The temporal evolution of chemical species, their lifetimes and relations between, for example, ozone and methane radiative forcing are analysed and compared to other studies in more detail in Sects.4.1 and 4.2.The climate impact metrics are derived by standard formulas (see e.g.Fuglestvedt et al., 2010) and a sanity check given in Sect.4.3.An overall evaluation of the metrics is given in Sect.4.3 by intercomparing the P-AGWP20 metrics derived from this work with results obtained with the climate-chemistry response model AirClim (Grewe and Stenke, 2008;Grewe and Dahlmann, 2012).

Verification
We have set up a complex modelling scheme, combining new methods for calculating the contribution of aviation on atmospheric processes with an air traffic model.A validation of this specific model application is not possible, since most of the effects are not yet measured or are per se not measurable, such as the effect of a local NO x emission at 200 hPa on the temperature change after 50 yr.Instead, we compare our results with earlier modelling studies.Even this is extremely difficult, since the focus of these studies is different and no direct intercomparison is possible.Hence we are more focusing on the soundness of the data (sanity check) rather than a verification or validation.An overview is given in Table 1.
The limited possibilities to directly compare the results of our modelling approach to either observational data or other model results emphasises the need to include sensitivity studies in future investigations on the changes in air traffic routing when optimising with respect to climate impact.This can be achieved by targeted manipulation of the climate cost functions with respect to, for example, the ratio of the impact from individual compounds, their variability or pattern and their consequences on the climate optimal routing changes.2), 50 air parcel trajectories are started and the temporal evolution of the global mean over these 50 air parcel trajectories is presented; that is, each coloured line (in total 504) represents the temporal evolution of species caused by an emission at the respective time-region grid point.The white lines indicate results from Stevenson et al. (2004) scaled by the factor 6.94 × 10 −3 to derive the same initial NO x perturbation.Note that values by Stevenson et al. (2004) are monthly means.

Chemistry
We have computed climate cost functions for one specific day in December.The underlying weather pattern represents a strong zonal jet, which represents winter pattern 1 according to Irvine et al. (2013).For each cost function grid point (see Table 2), the mean over all 50 air parcel trajectories is calculated.Figure 9 shows the temporal development of NO x , ozone and H 2 O as well as loss of CH 4 of these mean values for each time-region grid point.
A NO x and ozone lifetime of 20 ± 11 days and 72 ± 26 days, respectively, is calculated which is roughly in agreement with findings by Stevenson et al. (2004), who found a decrease of a NO x pulse emission from roughly 110 Gg(NO 2 ) to 20 Gg(NO 2 ) within a month, which represents a lifetime of approximately 25 days and a decrease of ozone from around 10 Tg to 1.5 Tg within 2 months, representing a lifetime of approximately 50 days (see also Fig. 9).The temporal evolution of NO x , ozone and methane changes for a January pulse emission calculated by Stevenson et al. (2004) (white lines in Fig. 9) is well within the range of our results.
The relation between the NO x emission and ozone contribution as a mean over all time-region grid points is 7.8 ± 2.4 DU(TgN a −1 ) −1 .Dahlmann et al. (2011), Frömming et al. (2012) and Grewe et al. (2002) calculated values for the whole air traffic and annual emissions of 1.4 DU(TgN a −1 ) −1 , 0.8 DU(TgN a −1 ) −1 and 0.7 DU(TgN a −1 ) −1 , respectively.However, those studies consider a whole air traffic scenario, with a large contribution of emissions from lower altitudes and lower ozone impacts.In our study, the ratio between the NO x enhancement and the ozone change is 2400 kg(O 3 ) per kg(N) with a range of 270 to 300 000 kg(O 3 ) per kg(N).Grewe et al. (2002) calculated a value of 300 kg(O 3 ) per kg(N); thus our result is again at the lower end for the same reason.The respective RF value is 41 mW m −2 DU −1 (this study) and compares well with 31 mW m −2 DU −1 calculated by Dahlmann et al. (2011).The emission-specific ozone RF calculated here ranges from 15 to 2800 mW m −2 TgN −1 with a mean value of 250 mW m −2 TgN −1 .Dahlmann et al. (2011) and Fuglestvedt et al. (2008) give mean values for the whole air traffic of 41 mW m −2 kgN −1 and 45 mW m −2 kgN −1 , respectively, which are well within the simulated range.
The relation between the ozone and methane RF is in the range of −0.5 to −1.3. Lee et al. (2010) and Holmes et al. (2011) summarised previous model simulations and found relations which amount to −1.65 ± 0.36 and −1.70, respectively.However, these values refer to the global air traffic, whereas we consider here emissions in the upper troposphere and lower stratosphere of the northern Atlantic region during winter.Grewe and Stenke (2008) give a value of around −1 at those altitudes and latitudes.
To summarise, the lifetime of simulated NO x and O 3 , the ratio between the NO x emission and ozone contribution, the resulting specific RF, and the ratio between the ozone and methane RF show large variabilities between the individual time-region grid points, and values published in the literature are well within this range.
The water vapour specific emission RF is 5.5 × 10 −11 mW m −2 (kg(H 2 O)a −1 ) −1 ranging from 0.5 to 20 × 10 −11 mW m −2 (kg(H 2 O)a −1 ) −1 .For a fleet with 170 Tg of fuel used per year this leads to an RF of 12 mW m −2 with a range of 1-43 mW m −2 .Here the value of Lee et al. (2010) of 3 mW m −2 is also well within this range.The recent estimate by Wilcox et al. (2012) of 0.9 mW m −2 with a range of 0.3 to 1.4 mW m −2 is at the lower end or slightly below of our extrapolation.Note again that their calculations were for the entire global fleet rather than for emissions in a particular latitude-height region, as is the case here.

Contrails
A variety of studies have investigated the properties of contrails (Lee et al., 2010;Heymsfield et al., 2011).Observations of contrails cover an age spectrum from seconds to hours and were performed by in situ measurement techniques and remote sensing, e.g. from satellite platforms.Here we focus on some climate-relevant contrail properties, such as ice water content and optical depth in the visible spectral range, and show a model-to-model radiative forcing benchmark test.
Both modelling and observational data vary by orders of magnitude (e.g.Kärcher et al., 2009;Schumann, 2012).In addition, direct intercomparisons are often challenging for many reasons, e.g.different environmental conditions, different sampling periods and different detection limits.For example, modelling studies often simulate a wide range of optical thickness of contrails, which, for example, cannot be detected from satellite (Marquart et al., 2003;Kärcher et al., 2009).For all these reasons, the comparison of contrail properties such as ice water content and optical thickness for a limited number of simulations can just be seen as a sanity check rather than a hard benchmark test (Grewe et al., 2012b).
Figure 10a shows observed and simulated ice water content in contrails as a cumulative probability density function.In-situ measurements from Voigt et al. (2011) (blue) and Schröder et al. (2000) (green) include 14 and 12 contrails, respectively, but nevertheless show a large variability.Our simulations (88 regions of contrail formation, red) cover a similar range and show a similar distribution to Voigt et al. (2011).Kärcher et al. (2009) simulated a wide range of environmental conditions for contrail formation and found a median ice water content of 0.6 mg m −3 in their modelling study, which is close to our results of 0.2 mg m −3 .
Note that the location and time of year differ in these studies.Data from Voigt et al. (2011) are obtained at 220 hPa to 300 hPa in central Europe in November and Schröder et al. (2000) include data in the altitude range from 300 hPa to 200 hPa also over central Europe, but in April/May and October, whereas here we investigate simulated contrail properties from 200 hPa to 400 hPa over the North Atlantic in December.We find that 88 out of 504 situations where a contrail forms, and take values of ice water content and optical properties when the contrail is fully evolved.In addition, the sampling frequency and representativity of the data with respect to the contrail area largely differs between the studies.
The cumulative probability density function of the contrail optical depth for wavelengths in the visible of these is presented in Fig. 10b.Satellite measurements (Iwabuchi et al., 2012) obligatorily deviate from in situ measurements, such as Voigt et al. (2011) and modelling studies such as Kärcher et al. (2009) and Frömming et al. (2011), because they cannot detect subvisible or thin contrails, not to mention that the sampling periods and areas are substantially different.In general, we find smaller optical depths for our specific simulated  We performed the RF-contrail benchmark test by Myhre et al. (2009), where a globally uniform 1 % contrail coverage with an optical depth of τ = 0.3 in the visible spectrum at 11 km is prescribed.The results are shown in Fig. 11 and 12. EMAC and ECHAM4 (E4) fall well within the variability of the all models (for details on this test for E4, see Frömming et al., 2011).EMAC has an increased number of spectral intervals compared to E4: 4 and 16 bands in the short-and longwave, compared to 2 and 6, respectively.The difference in the longwave part of the RF between E4 and EMAC is small, but it differs significantly for the shortwave net forcing.However, both models are well within the overall range, which is encouraging, as four of the five radiative transfer models of Myhre et al. (2009) are much more sophisticated than the simplified schemes that have to be used in 3-D models.

Metrics
Here we compare the results for three common climate metrics (AGWP20, AGWP100 and AGTP50) for ozone and methane (see Sect. 3.5) with results presented in Fuglestvedt et al. (2010).They analysed the GWP and GTP of ozone and methane (i.e. the AGWP and AGTP of ozone and methane normalised with the related AGWP and AGTP of CO 2 ) for different time horizons for a 1 yr pulse emission of NO x .We normalise the values with the respective metric for CO 2 (Table 8) to obtain the respective GWP and GTP values.For all metrics the ozone values are well within the range of previous studies.As already reported in Sect.4.1, the regarded region (mid-to high latitudes in the tropopause region) is characterised by a stronger methane loss per ozone enhancement than for the whole air traffic (Grewe and Stenke, 2008).Hence the methane metric values are at the lower end of the values reported in Fuglestvedt et al. (2010).
To obtain an overall assessment of the metric results (here P-AGWP20), we take the transatlantic air traffic emissions, calculated with the SAAM model for the minimum economic cost, and compare these results with an AirClim simulation based on the same emission data.AirClim is a fast climatechemistry response model (Grewe and Stenke, 2008;Grewe and Dahlmann, 2012) which takes into account annual mean emissions and their regional different effects (basically latitudes and altitudes) based on a number of precalculated cases  with complex chemistry-climate modelling.Note that in RE-ACT4C the climate response takes the specific weather situation into account, whereas in the case of the AirClim simulation, the emissions are assumed to occur everyday at the same place, i.e. identical for all weather situations throughout a year.The results (Fig. 13) show almost identical values for all emission components, except for contrail cirrus, which is reasonable, since the day-to-day variability of the contrail effects are highly variable and a one-day simulation cannot be expected to be representative of the whole year.However, the comparison shows that the overall results are comparable in magnitude.

Aircraft trajectory
As an illustration, Fig. 14 shows various flight options for one city pair connection (Washington to Vienna).The flight on that day (light brown) clearly follows the jet stream (arrows).
We have performed an optimisation of the whole transatlantic air traffic with respect to short-term climate impacts (question Q2 in Table 7, with option 3 GWP20) and costs, and we obtain a different aircraft trajectory for the above-mentioned city pair (blue), which more closely follows the jet stream.However, when the traffic is cost-optimised and includes conflict avoidance, as in reality, then the real route (light brown) and the cost-optimal (cheapest) route within conflictfree traffic (dark brown) are much closer.The difference in distance and time is around 1 % and in fuel about 3 %.For this city pair, the climate optimal routes (with and without conflict avoidance) with respect to short-term climate impact (see above) are located further north and at lower flight altitudes (FL330 and FL310).The conflict-free, climate-optimised route (green) avoids more contrails and leads to a decrease in contrail AGWP20 by 16 % and a decrease in NO x AGWP20 by 4 % with an increase in fuel consumption by 14 %, which is related to an increase in GWP20 from CO 2 by 1 %.This is one essentially random example for a particular weather pattern, a particular city pair and a particular choice of climate metrics.Beyond the fact that "climate-friendly" routes can indeed differ from the least-cost routes, no general conclusions can be drawn from this one example.It is the aim of REACT4C to produce a more systematic analysis across weather conditions and metrics.

Conclusions
We have developed a simulation framework for investigating climate change mitigation options for air traffic routing by avoiding climate-sensitive regions.It includes three major steps: the calculation of climate cost functions, the simulation and optimisation of air traffic according to these climate cost functions, and the estimation of the total mitigation gain.The climate cost functions describe the air traffic's contribution to climate change, which is caused by an emission at a certain time and location in the atmosphere.The processes we are regarding are ozone formation, methane loss, methane-induced ozone change, contrails (including the spread into cirrus), water vapour and carbon dioxide.
The simulation of physical and chemical processes is described using the chemistry-climate model EMAC, which we extended via two submodels AIRTRAC and CONTRAIL.They are described in detail in the Supplement (Frömming et al., 2013).By using the Lagrangian transport scheme ATTILA, we were able to simulate a climate cost function at a multitude of grid points in one simulation.We used state-of-the-art chemistry and microphysics for the simulation.New is the way in which the model is used and how the climate cost functions are calculated, which required a few new considerations, e.g.regarding the calculation of the radiative forcing.We developed a parameterisation, which relates the radiation changes caused by a 15 min pulse emission to an adjusted radiative forcing for ozone.We demonstrated the methodology for one specific simulated day in December which was characterised by a strong zonal jet stream over the Atlantic.We calculated the 4-D climate cost function data set and optimised the transatlantic air traffic with respect to economic costs and climate impact.Unfortunately a validation is hardly possible, since observational data for relations between emissions at flight altitude and climate impact is not available.Instead, we performed a sanity check of our calculated data and compared our results to previous studies, knowing that the comparison is very crude.In detail, we compared our calculated contribution of an emission at a very specific location and time to published results based on global and mostly annual air traffic emissions.Nevertheless, the ranges of, for example, contribution of NO x emissions to ozone and RF, contrail properties, etc., compare well with values from the literature.
It is important to note that uncertainties are associated with the calculation of the climate cost functions.For the optimisation, the absolute values of the climate costs are less important.More important is the relation of the climate impact of individual components and the spatial and temporal variability.In general, the uncertainty of a mean value and of any sensitivity is not necessarily correlated.Stevenson et al. (2006) showed in a multi-model intercomparison that the simulated ozone burden and methane lifetime have quite some variability but that the models very consistently simulate sensitivities.Grewe and Dahlmann (2012) intercompared results of the effect of flight altitude changes on ozone, water vapour and contrails based on Grewe and Stenke (2008), Köhler et al. (2008) and Rädel and Shine (2008) and found similar sensitivities for ozone and contrails, but larger ones for water vapour and methane.
The calculated climate cost functions will form the basis for a detailed analysis of the climate change mitigation potential of the air traffic system.First results indicate that a large potential already exists at present to significantly reduce the contribution of air traffic to anthropogenic climate change.In future publications we will investigate this potential in more detail.The optimisation with respect to climate might be an extreme scenario.Nonetheless, small changes to the air traffic system, i.e. the change in routing of a few very climate-sensitive flights, might already yield a large reduction of the air traffic's climate impact through reduction of contrail and NO x effects and with only small increases in fuel consumption.However, such conclusions are influenced greatly by the overall objective or, in other words, the aim of any climate change policy, which in turn controls the choice of the climate emission metric and the choice of time horizon for those metrics.Some choices would put a greater value on reducing the forcing due to short-lived components such as contrails and NO x , whilst others would put a greater value on reducing the forcing due to the longer-lived emissions, notably CO 2 .A necessity for the operational implementation of the methodology would be a reliable forecast of weather conditions in advance of the route-planning procedures and confidence that specific conditions, such as the occurrence of regions of ice supersaturation, could be forecast at the necessary level of detail.

Calculation of the adjusted radiative forcing for ozone
In this section we provide a methodology to calculate radiative forcings which will serve as an input to the climate metric calculations (Sect.3.5).Since we are considering relatively short pulse emissions, previous approaches have to be adapted.For all species, except for contrails, we consider the adjusted radiative forcing (e.g.Hansen et al., 1997).For ozone this requires new considerations, which are presented in the following sections.For contrails the difference between instantaneous and adjusted radiative forcing is small (Marquart, 2003) and neglected here.

A1 General approach for adjusted radiative forcing for ozone
Radiative forcing is the common metric to intercompare the global mean impacts of various components contributing to total climate change.As emphasised in previous works (Hansen et al., 1997;Stuber et al., 2001), the so-called stratosphere-adjusted radiative forcing (also known as the fixed dynamical heating approximation; Forster et al., 1997) is generally a better quantification than the instantaneous forcing for the climate impact that is to be expected from an ozone perturbation near or above the tropopause.Other definitions of radiative forcing are available (see e.g.Gregory et al., 2002), but these require lengthy climate model calculations and are not easily applicable to the relatively small radiative forcings due to aviation.The basic equation is which relates the global mean stratosphere-adjusted radiative forcing RF adj in [W m −2 ] linearly to the global mean equilibrium surface temperature response T eq surf in [K].RF adj is the radiative imbalance (at the tropopause or at the top of the atmosphere) induced by the forcing perturbation, determined after the stratosphere has re-adjusted to thermal equilibrium.λ in [K (W m −2 ) −1 ] is the climate sensitivity parameter which has been assumed, initially, to be a universal constant, independent of the nature or the distribution of the perturbation.If the instantaneous forcing RF inst instead of RF adj is used in Eq. (A1), λ becomes strongly perturbationdependent and even the sign of the radiative forcing may become inconsistent with the resulting surface temperature response in the case of an ozone perturbation at higher altitudes (e.g.around 20 km) (Hansen et al., 1997).The necessity to quantify perturbations in terms of RF adj poses a specific problem for the methodology used in the present study: RF adj cannot be directly calculated for a pulse perturbation, as such a perturbation is not stationary, meaning that the stratospheric temperature cannot adjust to a new equilibrium.In contrast, RF inst could easily be determined, even in such a case, at any location of the parcel in space and time, as no temperature adjustment is needed for calculating the instantaneous radiative flux change induced by a changing absorber.However, because we regard each of our short-lived perturbations as part of aviation climate impact as a whole (when stratospheric temperature adjustment are non-negligible), it has to be assessed by a metric consistent with RF adj .In the following we will discuss a thought experiment that illustrates the problem.The general idea is to translate an instantaneous radiative forcing RF inst resulting from a pulse emission into an equivalent stratosphere-adjusted radiative forcing RF adj by means of an analytical formula: The unitless functions f 1 and f 2 describe the relation between RF adj and RF inst for various times of the year (f 1 , seasonal cycle) and various perturbation altitudes (f 2 ) and are derived from additional idealised simulations.

A2 Thought experiment
Let us assume we consider two alternative flight routings, implying emissions in two different time regions.We want to answer the following question: "which routing induces the lower climate impact?"Let us further assume we had the same meteorology everyday.The constraint from this thought experiment is then that the decision on the preferable routing is identical everyday.
We now investigate two different assessment approaches, which both should lead to the same choice of routing.Figure A1 shows RF inst (green) of a short pulse emission, e.g. as originating from NO x emissions of one individual flight.The pulse is short compared to the temperature adjustment time of the stratosphere and a temperature adjustment is not yet achieved.One can regard this period as a spin-up for RF adj , during which the adjusted and instantaneous RFs are still equal.Both RF inst and RF adj are positive in this case.A sequence of individual pulses (flights), occurring every 0.005 time units, will also lead to a positive, though higher, RF inst (blue).In this case, however, the continuous sequence of pulses gradually induces stratospheric temperature adjustment, reducing RF adj with respect to RF inst over time until it eventually becomes -in this example, which is geared to ozone effects -negative (−2 RF units, magenta).Thus, the assessment of an individual flight as part of a sequence of the same flights, and with the same meteorology everyday, has to  A1.Schematic view of the radiative forcing due to ozone as a result of a pulse emission of NO x compared to a sequence of pulse emissions (i.e.sustained emissions).Instantaneous and adjusted RF of a short pulse emission are identical (green).Instantaneous RF of a sequence of pulse emissions (blue) is calculated as the sum of the instantaneous RF of pulse emissions.Adjusted RF of a sequence of pulse emissions (magenta) is calculated similarly, i.e. as the sum of the RF of a sequence of pulse emissions, where the cooling effect of preceding pulses is parameterised (red).reflect this negative RF.Its time developing radiative impact consists of contributions of varying strength and sign, and it contains an additional tail of negative values (red line).This contribution must be included for the assessment of the flight in terms of its time-averaged RF.
As mentioned, such an equivalent RF adj cannot be gained directly from the individual pulse simulation.We overcome this dilemma by calculating a response function (Eq.A2) to obtain a parameterised RF adj , transferring RF inst to RF adj by means of precalculated temperature adjustment effects from various types of pulse perturbation.The terms of the response functions can be obtained from idealised simulations that calculate RF adj and the associated stratospheric temperature adjustment from a limited number of typical ozone perturbations located at various altitudes and latitudes, and under various conditions of insolation.

A3 Simulation setup
Three idealised ozone perturbation patterns (see Supplement) are considered, which were adopted from previous simulations evaluating a number of emissions, located in various latitude-height bands (Fichter, 2009).The NO x emissions basic to the ozone patterns discussed here occurred at altitudes of 200, 160 and 130 hPa, within the Northern Hemisphere extratropics (see Supplement for more details).Radiative forcing calculations were performed with EMAC for each perturbation pattern over a 1 yr period, preceded by a 3-month spin-up for the adjustment of stratospheric temperatures.The simulations include a calculation of the adjusted as well as the instantaneous RF.Each of the 12 months is interpreted as an individual pulse, and the RF inst and RF adj that it induces can be directly compared from the simulation, and the difference between RF inst and RF adj can be determined for each calender month.
Additionally, we have extended the simulation by four simulations for the second year, where we have successively switched off the ozone perturbation after January, April, July and October.These simulations are used to quantify the contribution of the adjusted stratospheric temperatures to the RF calculations after the perturbation has faded out.As soon as the perturbation pattern is removed in the respective (second simulation) month, the only contribution to the RF arises from the remaining stratospheric temperature changes, which revert to the unperturbed situation rather fast.It turned out that this contribution to RF adj is only of minor importance and thus we could omit it for the sake of simplicity.

A4 Separation of height and time dependencies
Figure A2a shows the seasonal cycle of the radiative forcing (both RF inst and RF adj ) as caused by a NO x pulse emission over one month (see Supplement for the respective ozone pattern).As explained in Sect.A3, RF adj values imply the assumption that preceding emissions induced a stratospheric temperature adjustment before the considered pulse was initiated.A distinct annual cycle is obvious, reflecting the insolation variation that directly effects the shortwave component  Stuber (2003) and Fichter (2009).Ozone changes are mainly found in the Northern Hemisphere (see also the Supplement).
of the radiative forcing.The contribution of the stratospheric temperature adjustment to RF adj also displays a seasonal cycle, as heating due to solar absorption by stratospheric ozone forms an important component of the adjustment.There is a general increase of RF with the altitude of the pulse emission.
For pulses emitted at 130 hPa and 160 hPa, RF adj is larger than RF inst , whereas the opposite is true for an emission at 200 hPa.Following our approach to derive a scaled instantaneous forcing serving as a proxy for RF adj , we normalise all seasonal cycles to the respective August value (Fig. A2b).The ratio between these August reference values of RF adj and RF inst is 1.56, 1.18 and 0.97 for 130 hPa, 160 hPa and 200 hPa, respectively.The relative time development becomes much more unique, and the normalised, adjusted RF is always lower than the normalised, instantaneous RF, allowing for the rescaling of RF inst to yield the required proxy RF adj (Eq.A2).The relative difference displays a clear seasonal cycle (Fig. A2c) with largest differences in winter (about 15 %), which can be approximated by a sine function (black line): with f mean 1 = 0.08 being the amplitude describing the deviation of peak values at summer and winter from the mean and t being the month of the year.It is obvious from Fig. A2c that the approximation is largely independent of the emission altitude, justifying the approach of decomposed functions in Eq. (A2).
The altitude dependency of the relation between RF inst and RF adj can be described from the three examples in Fig. A2a (based on simulations by Fichter, 2009, as mentioned)  evaluated by Stuber (2003).The fifth data point is at 50 hPa, and hence out of the region of interest, but used for the fit forming the parameterisation.Note that the experimental designs of Fichter (2009) and Stuber (2003) are both Northern Hemisphere ozone changes at different altitudes, though they differ in detail, especially for the prescribed ozone pattern.The results are displayed in Fig. A3.In cases where the ozone increase is exclusively located in the troposphere (i.e.below approximately 200 hPa), RF adj is smaller than RF inst , since the stratospheric temperature adjustment implies a cooling caused by the blanketing effect of the perturbation.In the tropopause region and in the lowermost stratosphere, an ozone increase induces a dipole-like stratospheric temperature adjustment, with warming by solar absorption where the stratospheric part of the perturbation peaks and cooling above that.The cases are also characterised by a positive instantaneous forcing at the tropopause, with longwave and shortwave contributions adding constructively (for tropospheric ozone changes), or with the positive longwave contribution dominating over a negative shortwave contribution (lowermost stratosphere).In cases where the ozone increase takes place at even higher altitudes (between 150 hPa and 50 hPa), RF inst at the tropopause becomes negative, because the effect of shortwave absorption dominates the instantaneous longwave cooling (Hansen et al., 1997;Stuber et al., 2001).RF adj continues to be positive, however, as the warming of the stratosphere by the absorption of the upwelling longwave radiation and downwelling shortwave radiation provides an additional downward longwave flux, which for lower tropospheric ozone perturbations is strong enough to overcompensate for the instantaneous net effect.In those cases the ratio between RF adj and RF inst is negative (Stuber, 2003).The fit resembles the overall structure well: with A = 1.1, B = 2.0, C = 1.05,D = 3.4,E = 1/60.A summarising sketch is provided in Fig. A4: the instantaneous upward and downward net flux (net means sum of shortwave and longwave) are balanced (black arrows) in the unperturbed case (left).When introducing a tropospheric perturbation (second left, red line) the downward fluxes remain basically unchanged, but the upward longwave flux (and hence the upward net flux) is reduced.This instantaneous flux imbalance (green arrow) equals a positive instantaneous radiative forcing (green bar).This leads to a stratospheric cooling (blue line) and hence to a reduced downward longwave flux, shown as an upward net flux (blue arrow).The combination of the instantaneous radiative forcing and this flux change yields the adjusted radiative forcing, which is hence smaller than the instantaneous RF.
A perturbation near tropopause altitudes (second right, red line) is much more effective than a pure tropospheric perturbation, leading to a stronger flux imbalance (green arrow) and instantaneous RF (green bar).The adjusted stratospheric temperature, however, has a cooling and warming component (blue line) and may result in an additional longwave downward flux (blue arrow) and hence a larger adjusted RF than instantaneous RF.A pure stratospheric perturbation (right) leads to a reduced shortwave downward flux and hence a reduced net flux (black arrow).The instantaneous RF is therefore negative, while the resulting adjusted stratospheric temperatures are positive (blue line).This leads to a strong longwave downward flux (blue arrow), which overcompensates for the flux imbalance (green arrow) and yields a positive adjusted RF.

Fig. 6 .
Fig.6.Top: deviations [%] of possible results of the monthly mean January NO x mixing ratio (left) and ozone mixing ratio (right) for emissions at 200 hPa and 50 • N, 45 • W as a function of the number of air parcel trajectories simulated.The deviation is relative to the result using 50 air parcel trajectories.Bottom: standard deviations [%] of the results for NO x and ozone shown as a function of the number of air parcel trajectories used for 24 different time-region grid points (shading: grey for NO x and blue for ozone) and the mean of the standard deviation over the 24 points (upper black line for NO x , lower dark-blue line for ozone).
In an analogous way, we derive the production and loss terms for HO e 2 and obtain differential equations for OH e and HO e 2 , which easily can be solved as follows: www.geosci-model-dev.net/7/175/2014/Geosci.Model Dev., 7, 175-201, 2014 V. Grewe et al.: Modelling of climate cost functions

Fig. 8 .
Fig. 8. Relation between change in atmospheric water vapour [Tg] and adjusted radiative forcing [mW m −2 ] for a fixed emission at different locations.Note that this relation is only valid if a constant emission is taken at different locations, which results in both different mass changes and different RFs.Data are from Grewe and Stenke (2008).The line represents a fit with RF = 4.38 × 10 −13 W m −2 kg −1 .

Fig. 9 .
Fig. 9. Temporal evolution of global mean masses of NO x [10 3 kg(NO)], O 3 [10 6 kg], H 2 O [10 5 kg] and CH 4 [10 5 kg].For every timeregion grid point (see Table2), 50 air parcel trajectories are started and the temporal evolution of the global mean over these 50 air parcel trajectories is presented; that is, each coloured line (in total 504) represents the temporal evolution of species caused by an emission at the respective time-region grid point.The white lines indicate results fromStevenson et al. (2004) scaled by the factor 6.94 × 10 −3 to derive the same initial NO x perturbation.Note that values byStevenson et al. (2004) are monthly means.

Fig. 10 .
Fig. 10.Intercomparison of observed and simulated distributions of contrail properties.Estimated cumulative probability density functions of contrail ice water content [mg m −3 ] (a) and (b) visible optical depth.

Fig. 11 .
Fig. 11.Geographical distribution of the annual mean all-sky net radiative forcing at the top of the atmosphere for a homogeneous 1 % contrail cover simulated with EMAC according to the benchmark test by Myhre et al. (2009).
day than the studies byKärcher et al. (2009) andFrömming et al. (2011), which cover a much broader range of environmental parameters.Kärcher et al. (2010) compared simulated optical depths with a cloud model for the United States and found a median visible optical depth of 0.02, which is closer to our simulated median visible optical depth.

Fig. 12 .
Fig. 12. Global means for the longwave (top), shortwave (mid) and net RF (bottom) for the contrail RF benchmark test (see text for details).Figure is adapted from Myhre et al. (2009).

Fig. 13 .
Fig. 13.Comparison of the pulse absolute global warming potential (P-AGWP20) in [µW/m 2 ] for emissions of CO 2 (red), contrail cirrus (green), NO x (blue) and H 2 O (magenta).Note that the H 2 O values are multiplied by 100 for presentational purposes.The emissions from the minimum economic cost one-day transatlantic air traffic, which are calculated with SAAM, are multiplied with the climate cost functions (REACT4C, left) and taken as annual mean emissions for the climate-chemistry response model AirClim (right).

Fig. 14 .
Fig. 14.Examples of climate-optimised flights.The optimisation was performed with respect to question Q1 (short-term climate impacts, AGWP20, see Table7) for one particular winter weather pattern (zonally strong jet).The selected city pair connection is Washington to Vienna.The real flight at that day is shown in light brown.The economic optimal flights without conflict avoidance is given in blue and the conflict avoidance is given in dark brown.The climate optimal flights are shown in green.In this case the aircraft trajectories with and without conflict avoidance differ only in cruise altitude.Arrows indicate the wind field at flight level 380, i.e. 38 000 ft.

Fig. A2 .
Fig. A2.Annual cycle of the radiative forcing due to ozone change of a pulse emission of NO x .The NO x emission which produces the monthly change pattern, is located at 130 hPa (L11, red), 160 hPa (L13, blue) and 200 hPa (L15, green).Dashed lines refer to instantaneous RF and solid ones to adjusted RF.(a) Seasonal cycle of RF.(b) as (a), but normalised to the individual August value.(c) Relative difference [%] of the adjusted RF and the scaled instantaneous RF.The scaling factor is the quotient of the respective August values.The black line shows a sine fit f 1 (see Eq. (A3)).
Fig. A3.Height dependence of the relation between instantaneous and adjusted RF for annual mean RF values of changes in ozone.Data are taken fromStuber (2003) andFichter (2009).Ozone changes are mainly found in the Northern Hemisphere (see also the Supplement).
Fig. A4.Schematic of the relation between instantaneous (green bar) and adjusted RF (magenta bar) for ozone changes for four situations: unperturbed, tropospheric perturbation, tropopause perturbation and stratospheric perturbation (from left to right).The perturbation pattern is given as a red line, the instantaneous upward and downward net flux changes are given as black arrows.Net flux means the sum of longwave and shortwave fluxes.The adjusted stratospheric temperatures are shown as a blue line, and the resulting net flux changes as blue arrows.For details characterising the different cases, see the text.

Table 1 .
Overview on the climate cost function calculation, assumptions and verification.

Table 2 .
Definition of the time-region grid.

Table 3 .
Emissions and aircraft/fuel parameters.

Table 4 .
Parameters of the contrail parameterisation.

Table 5 .
Parameters of the CO 2 parameterisation.

Table 6 .
Overview on radiative forcing calculation for various species."A" and "I" stand for adjusted and instantaneous radiative forcing in the column "RF type".

Table 7 .
Overview of objectives and implied emission scenarios and climate metrics or indicators."F","S" and "P" indicate assumptions on emissions: future emission scenario, sustained emissions and pulse emissions, respectively.ATR is the average temperature response, AGTP the absolute global temperature potential, and AGWP the absolute global warming potential; the number after the metrics gives the time horizon in years.Q1: What is the short-term climate impact Q2: What is the long-term climate impact Q3: What is the medium-range climate impact

Table 8 .
Fuglestvedt et al. (2010)WP20, GWP100 and GTP50 values for ozone (short-lived only) and methane for emissions of NO x in the domain indicated in Table2with results from three different studies published inFuglestvedt et al. (2010), abbreviated as F10.