Air traffic simulation in chemistry-climate model EMAC 2.41

. Mobility is becoming more and more important to society and hence air transportation is expected to grow further over the next decades. Reducing anthropogenic climate impact from aviation emissions and building a climate-friendly air transportation system are required for a sustainable development of commercial aviation. A climate optimized routing, which avoids climate-sensitive regions by re-routing horizontally and vertically, is an important measure for climate impact reduction. The idea includes a number of different routing strategies (routing options) and shows a great potential for the reduction. To evaluate this, the impact of not only CO 2 but also non-CO 2 emissions must

Abstract.Mobility is becoming more and more important to society and hence air transportation is expected to grow further over the next decades.Reducing anthropogenic climate impact from aviation emissions and building a climatefriendly air transportation system are required for a sustainable development of commercial aviation.A climate optimized routing, which avoids climate-sensitive regions by rerouting horizontally and vertically, is an important measure for climate impact reduction.The idea includes a number of different routing strategies (routing options) and shows a great potential for the reduction.To evaluate this, the impact of not only CO 2 but also non-CO 2 emissions must be considered.CO 2 is a long-lived gas, while non-CO 2 emissions are short-lived and are inhomogeneously distributed.This study introduces AirTraf (version 1.0) that performs global air traffic simulations, including effects of local weather conditions on the emissions.AirTraf was developed as a new submodel of the ECHAM5/MESSy Atmospheric Chemistry (EMAC) model.Air traffic information comprises Eurocontrol's Base of Aircraft Data (BADA Revision 3.9) and International Civil Aviation Organization (ICAO) engine performance data.Fuel use and emissions are calculated by the total energy model based on the BADA methodology and Deutsches Zentrum für Luft-und Raumfahrt (DLR) fuel flow method.The flight trajectory optimization is performed by a genetic algorithm (GA) with respect to a selected routing option.In the model development phase, benchmark tests were performed for the great circle and flight time routing options.
The first test showed that the great circle calculations were accurate to −0.004 %, compared to those calculated by the Movable Type script.The second test showed that the optimal solution found by the algorithm sufficiently converged to the theoretical true-optimal solution.The difference in flight time between the two solutions is less than 0.01 %.The dependence of the optimal solutions on the initial set of solutions (called population) was analyzed and the influence was small (around 0.01 %).The trade-off between the accuracy of GA optimizations and computational costs is clarified and the appropriate population and generation (one iteration of GA) sizing is discussed.The results showed that a large reduction in the number of function evaluations of around 90 % can be achieved with only a small decrease in the accuracy of less than 0.1 %.Finally, AirTraf simulations are demonstrated with the great circle and the flight time routing options for a typical winter day.The 103 trans-Atlantic flight plans were used, assuming an Airbus A330-301 aircraft.The results confirmed that AirTraf simulates the air traffic properly for the two routing options.In addition, the GA successfully found the time-optimal flight trajectories for the 103 airport pairs, taking local weather conditions into account.The consistency check for the AirTraf simulations confirmed that calculated flight time, fuel consumption, NO x emission index and aircraft weights show good agreement with reference data.
Published by Copernicus Publications on behalf of the European Geosciences Union.

Introduction
World air traffic has grown significantly over the past 20 years.With the increasing number of aircraft, the air traffic's contribution to climate change becomes an important problem.Nowadays, aircraft emission (this includes still uncertain aviation-induced cirrus cloud effects) contributes approximately 4.9 % (with a range of 2-14 %, which is a 90 % likelihood range) of the total anthropogenic radiative forcing (Lee et al., 2009;Lee et al., 2010;Burkhardt and Kärcher, 2011).An Airbus forecast shows that the world air traffic might grow at an average annual rate of 4.6 % over the next 20 years (2015( -2034( , Airbus, 2015)), while Boeing forecasts a value of 4.9 % over the same period (Boeing, 2015).This implies a further increase in aircraft emissions and therefore environmental impacts from aviation rise.Reducing the impacts and building a climate-friendly air transportation system are required for a sustainable development of commercial aviation.The emissions induced by air traffic primarily comprise carbon dioxide (CO 2 ), nitrogen oxides (NO x ), water vapor (H 2 O), carbon monoxide, unburned hydrocarbons and soot.They lead to changes in the atmospheric composition, thereby changing the greenhouse gas concentrations of CO 2 , ozone (O 3 ), H 2 O and methane (CH 4 ).The emissions also induce cloudiness via the formation of contrails, contrail-cirrus and soot cirrus (Penner et al., 1999).
The climate impact induced by aircraft emissions depends partially on local weather conditions.That is, the impact depends on geographical location (latitude and longitude) and altitude at which the emissions are released (except for CO 2 ) and time.In addition, the impact on the atmospheric composition has different timescales: chemical effects induced by the aircraft emissions have a range of lifetimes and affect the atmosphere from minutes to centuries.CO 2 has long perturbation lifetimes on the order of decades to centuries.The atmosphere-ocean system responds to the change in the radiation fluxes on the order of 30 years.NO x , released in the upper troposphere and lower stratosphere, has a different lifetime ranging from a few days to several weeks, depending on atmospheric transport and chemical background conditions.In some regions, which experience a downward motion, e.g., ahead of a high-pressure system, NO x has short lifetimes and is converted to HNO 3 and then rapidly washed out (Matthes et al., 2012;Grewe et al., 2014b).The most localized and short-lived effect is contrail formation, with typical lifetimes from minutes to hours.Persistent contrails only form in ice supersaturated regions (Schumann, 1996) and extend a few hundred meters vertically and about 150 km along a flight path (with a standard deviation of 250 km) with a large spatial and temporal variability (Gierens and Spichtinger, 2000;Spichtinger et al., 2003).
The measures to counteract the climate impact induced by aircraft emissions can be classified into two categories: technological and operational measures, as summarized by Irvine et al. (2013).The former includes aerodynamic im-provements of aircraft (blended-wing-body aircraft, laminar flow control, etc.), more efficient engines and alternative fuels (liquid hydrogen, bio-fuels).The latter includes efficient air traffic control (reduced holding time, more direct flights, etc.), efficient flight profiles (continuous descent approach) and climate-optimized routing.Nowadays, flight trajectories are optimized with respect to time and economic costs (fuel, crew, other operating costs) primarily by taking advantage of tail winds, e.g., jet streams, while the climate-optimized routing should optimize flight trajectories such that released aircraft emissions lead to a minimum climate impact.Earlier studies investigated the effect of systematic flight altitude changes on the climate impact (Koch et al., 2011;Schumann et al., 2011;Frömming et al., 2012;Søvde et al., 2014).They confirmed that the changed altitude has a strong effect on the reduction of climate impact.A number of studies have investigated the potential of applying climate-optimized routing for real flight data.Matthes et al. (2012) and Sridhar et al. (2013) addressed weather-dependent trajectory optimization using real flight routes and showed a large potential of climate-optimized routing.As the climate impact of aircraft emissions depends on local weather conditions, Grewe et al. (2014a) optimized flight trajectories by considering regions described as climate-sensitive regions and showed a trade-off between climate impact and economic costs.That study reported that large reductions in the climate impact of up to 25 % can be achieved by only a small increase in economic costs of less than 0.5 %.The climate-optimized routing therefore seems to be an effective routing option for the climate impact reduction; however, this option is still unused in today's flight planning.
This paper presents the new AirTraf submodel (version 1.0, Yamashita et al., 2015) that performs global air traffic simulations coupled to the EMAC chemistry-climate model (Jöckel et al., 2010).This paper technically describes Air-Traf and validates the various components for simple aircraft routings: great circle and time-optimal routings.Eventually, we are aiming at an optimal routing for climate impact reduction.The development described in this paper is a prerequisite for the investigation of climate-optimized routings.The research road map for our study is as follows (Grewe et al., 2014b): the first step was to investigate the influence of specific weather situations on the climate impact of aircraft emissions (Matthes et al., 2012;Grewe et al., 2014b).This results in climate cost functions (CCFs, Frömming et al., 2013;Grewe et al., 2014a, b) that identify climate sensitive regions with respect to O 3 , CH 4 , H 2 O and contrails.They are specific climate metrics, i.e., climate impacts per unit amount of emission, and will be used for optimal aircraft routings.In a further step, weather proxies will be identified for the specific weather situations, which correlate the intensity of the climate-sensitive regions with meteorological data.The proxies will be available from numerical weather forecasts, like temperature, precipitation, ice-supersaturated regions, vertical motions or weather patterns in general.These prox-ies are then used to optimize air traffic with respect to the climate impact expressed by the CCFs.An assessment platform is required to validate the optimization strategy based on the proxies in multi-annual (long-term) simulations and to evaluate the total mitigation gain of the climate impactone important objective of the AirTraf development.
This paper is organized as follows.Section 2 presents the model description and calculation procedures of AirTraf.Section 3 describes aircraft routing methodologies for the great circle and flight time routing options.A benchmark test provides a comparison of resulting great circle distances with those calculated by the Movable Type script (MTS, Movable Type script, 2014).Another benchmark test compares the optimal solution to the true-optimal solution.The dependence of optimal solutions on initial populations (a technical terminology set in italics is explained in Table A1 in Appendix A) is examined and the appropriate population and generation sizing is discussed.In Sect.4, AirTraf simulations are demonstrated with the two options for a typical winter day (called 1-day AirTraf simulations) and the results are discussed.Section 5 verifies whether the AirTraf simulations are consistent with reference data and Sect.6 concludes the study.Finally, Sect.7 describes the code availability.
2 AirTraf: air traffic in a climate model

Overview
AirTraf was developed as a submodel of EMAC (Jöckel et al., 2010) to eventually assess routing options with respect to climate.This requires a framework where we can optimize routings every day and assess them with respect to climate changes.EMAC provides an ideal framework, since it includes various submodels, which actually evaluate climate impact, and it simulates local weather situations on long timescales.As stated above, we were focusing on the development of this model.A publication on the climate assessment of routing changes will be published as well.
Figure 1 shows the flowchart of the AirTraf submodel.First, air traffic data and AirTraf parameters are read in messy_initialize, which is one of the main entry points of the Modular Earth Submodel System (MESSy, Fig. 1, dark blue).Second, all entries are distributed in parallel following a distributed memory approach (messy_init_memory, Fig. 1, blue): AirTraf is parallelized using the message passing interface (MPI) standard.As shown in Fig. 2, the 1-day flight plan, which includes many flight schedules of a single day, is decomposed for a number of processing elements (PEs; here PE is synonymous with MPI task), so that each PE has a similar workload.A whole flight trajectory between an airport pair is handled by the same PE.Third, a global air traffic simulation (AirTraf integration, Fig. 1, light blue) is performed in messy_global_end, i.e., at the end of the time loop of EMAC.Thus, both short-term and long-term simulations can take into account the local weather conditions for every flight.This AirTraf integration is linked to several modules: the aircraft routing module (Fig. 1, light green) and the fuel/emissions calculation module (Fig. 1, light orange).The former is also linked to the flight trajectory optimization module (Fig. 1, dark green) to calculate flight trajectories corresponding to a selected routing option.The latter calculates fuel use and emissions on the calculated trajectories.Finally, the calculated flight trajectories and global fields (three-dimensional emission fields) are output (Fig. 1, rose red).The results are gathered from all PEs for output.The output will be used to evaluate the reduction potential of the routing option on the climate impact.
The following assumptions are made in AirTraf (version 1.0): a spherical Earth is assumed (radius is R E = 6371 km).The aircraft performance model of Eurocontrol's Base of Aircraft Data (BADA Revision 3.9, Eurocontrol, 2011) is used with a constant Mach number M (the Mach number is the velocity divided by the speed of sound).When an aircraft flies at a constant Mach number, the true air speed V TAS and ground speed V ground vary along flight trajectories.Only the cruise flight phase is considered, while ground operations, take-off, landing and any other flight phases are unconsidered.Potential conflicts of flight trajectories and operational constraints from air traffic control, such as the semi-circular rule (the basic rule for flight level) and limit rates of aircraft climb and descent, are disregarded.However, a workload analysis of air traffic controllers can be performed on the basis of the output data.The following sections describe the used models briefly, while characteristic procedures of AirTraf are described in detail.

Chemistry-climate model EMAC
The ECHAM/MESSy Atmospheric Chemistry (EMAC) model is a numerical chemistry and climate simulation system that includes submodels describing tropospheric and middle atmosphere processes and their interaction with oceans, land and influences coming from anthropogenic emissions (Jöckel et al., 2010).It uses the second version of the MESSy (i.e., MESSy2) to link multi-institutional computer codes.The core atmospheric model is the 5th generation European Centre Hamburg general circulation model (ECHAM5, Roeckner et al., 2006).For the present study we applied EMAC (ECHAM5 version 5.3.02,MESSy version 2.41) in the T42L31ECMWF resolution, i.e., with a spherical truncation of T42 (corresponding to a quadratic Gaussian grid of approximately 2.8 • by 2.8 • in latitude and longitude) with 31 vertical hybrid pressure levels up to 10 hPa (middle of the uppermost layer).MESSy provides interfaces (Fig. 1, yellow) to couple various submodels.Further information about MESSy, including the EMAC model system, is available from http://www.messy-interface.org.  .Flowchart of EMAC/AirTraf.MESSy as part of EMAC provides interfaces (yellow) to couple various submodels for data exchange, run control and data input/output.Air traffic data and AirTraf parameters are input in the initialization phase (messy_initialize, dark blue).AirTraf includes the flying process in messy_global_end (dashed box, light blue), which comprises four main computation procedures (bold-black boxes).The detailed procedures are described in Sect.2.4 and are illustrated in Fig. 3. AirTraf is linked to three modules: the aircraft routing module (light green), the flight trajectory optimization module (dark green), and the fuel/emissions calculation module (light orange).Resulting flight trajectories and global fields are calculated for output (rose red).Various submodels of EMAC can be linked to evaluate climate impacts on the basis of the output.A whole trajectory of an airport pair is handled by the same PE in the time loop of EMAC (messy_global_end, light blue).Finally, results are gathered from all the PEs for output (rose red).

Air traffic data
The air traffic data (Fig. 1, dark blue) consist of a 1-day flight plan and aircraft and engine performance data.Concerning the engine performance data, four data pairs of reference fuel flow f ref (in kg(fuel) s −1 ) and the corresponding NO x emission index EINO x,ref (in g(NO x ) (kg(fuel)) −1 ) at take-off, climb out, approach and idle conditions are taken from the ICAO engine emissions databank (ICAO, 2005).An overall (passenger/freight/mail) weight load factor is also provided by ICAO (Anthony, 2009).

Calculation procedures of the AirTraf submodel
The calculation procedures in the AirTraf integration are described here step by step.As shown in Fig. 1 (light blue), the flight status of all flights is initialized as "non-flight" at the first time step of EMAC.The departure check is then performed at the beginning of every time step.When a flight gets to the time for departure in the time loop of EMAC, its flight status changes into "in-flight".The time step index of EMAC t is introduced here.The index is assigned t = 1 to the flight at the departure time.Thereafter the flight moves to the flying process (dashed box in Fig. 1, light blue), which mainly comprises four steps (bold-black boxes in Fig. 1, light blue): flight trajectory calculation, fuel/emissions calculation, aircraft position calculation and gathering global emissions.The following parts of this section describe these four steps and Fig. 3a to d illustrate the respective steps.
The flight trajectory calculation linked to the aircraft routing module (Fig. 1, light green) calculates a flight trajectory corresponding to a routing option.AirTraf will provide seven routing options: great circle (minimum flight distance), flight time (time-optimal), NO x , H 2 O, fuel (which might differ from H 2 O, if alternative fuel options can be used), contrail and CCFs (Frömming et al., 2013;Grewe et al., 2014b).In AirTraf (version 1.0), the great circle and the flight time routing options can currently be used.The great circle option is a basis for the other routing options and the module calculates a great circle by analytical formulae, assuming constant flight altitude.In contrast to this, for the other six op-tions, a single-objective minimization problem is solved for the selected option by the linked flight trajectory optimization module (Fig. 1, dark green); this module comprises the Genetic Algorithm (GA, Holand, 1975;Goldberg, 1989) and finds an optimal flight trajectory including altitude changes.For example, if the flight time routing option is selected, the flight trajectory optimization is applied to all flights taking into account the individual departure times.Generally, a wind-optimal route means an economically optimal flight route taking the most advantageous wind pattern into account.This route minimizes total costs with respect to time, fuel and other economic costs; i.e., it has multiple objectives.AirTraf will provide the flight time and the fuel routing options to investigate trade-offs (conflicting scenarios) among different routing options.With the contrail option, the best trajectory for contrail avoidance will be found.The CCFs are provided by EU FP7 project REACT4C (Reducing Emissions from Aviation by Changing Trajectories for the benefit of Climate, REACT4C, 2014) and estimate climate impacts due to some aviation emissions (see Sect. 1).Thus, the best trajectory for minimum CCFs will be calculated.
For all routing options, local weather conditions provided by EMAC at t = 1 (i.e., at the departure day and time of the aircraft) are used to calculate the flight trajectory.The con-  ditions are assumed to be constant during the flight trajectory calculation.No weather forecasts (or weather archives) are used.Once an optimal flight trajectory is calculated, it is not re-optimized in subsequent time steps (t ≥ 2).The detailed flight trajectory calculation methodologies for the great circle and the flight time routing options are described in Sect.3.After the flight trajectory calculation, the trajectory consists of waypoints generated along the trajectory, and flight segments (Fig. 3a).In addition, a number of flight properties are available corresponding to the waypoints, flight segments and the whole trajectory, as listed in Table 2. Here, the waypoint index i is introduced (i = 1, 2, • • •, n wp ); n wp is the number of waypoints arranged from the departure airport (i = 1) to the arrival airport (i = n wp ).i is also used as the flight segment index Next, fuel use, NO x and H 2 O emissions are calculated by the dedicated module (Fig. 1, light orange); this module comprises a total energy model based on the BADA methodology (Schaefer, 2012) and the DLR fuel flow method (Deidewig et al., 1996; see Sects.2.5 and 2.6 for more details).After this calculation, additional flight properties are newly available (see Fig. 3b and Table 2).Note that the flight trajectory calculation described above and this fuel/emissions calculation are performed only once at t = 1.
Table 2. Properties assigned to a flight trajectory.The properties of the three groups (divided by rows) are obtained from the nearest grid point of EMAC, the flight trajectory calculation (Fig. 3a), and the fuel/emissions calculation (Fig. 3b), respectively.The attribute type indicates where the values of properties are allocated."W", "S" and "T" stand for waypoints and a whole flight trajectory in column 3, respectively.

Property Unit
Attribute type Description The next step is to advance the aircraft positions along the flight trajectory corresponding to the time steps of EMAC (Fig. 3c).Here, aircraft position parameters pos new and pos old are introduced to indicate the present position (at the end of the time step) and previous position (at the beginning of the time step) of the aircraft along the flight trajectory.They are expressed by real numbers as a function of the waypoint index i (integers), i.e., real(1, 2, • • •, n wp ).At t = 1, the aircraft is set at the first waypoint (pos new = pos old = 1.0).As the time loop of EMAC progresses, the aircraft moves along the trajectory referring to the estimated time over (ETO, Table 2) (AirTraf continuously treats overnight flights with arrival on the next day).For example, Fig. 3c shows pos new = 2.3 and pos old = 1.0 at t = 2.This means that the aircraft moves 100 % of the distance between i = 1 and i = 2, and 30 % of the distance between i = 2 and i = 3 in one time step.pos new and pos old are stored in the memory and the aircraft continues the flight from pos new = 2.3 at the next time step.After the aircraft moves to a new position, the arrival check is performed (dashed box in Fig. 1, light blue).If pos new ≥ real(n wp ), the flight status changes to "arrived".
Finally, the individual aircraft's emissions corresponding to the flight path in one time step are gathered into a global field (three-dimensional Gaussian grid).This step is applied for all flights with "in-flight" or "arrived" status.As shown in Fig. 3d, for example, the released NO x emission along a flight segment i (NO x,i or the fraction of it) is mapped onto the nearest grid point of the global field.For this NO x,i , the coordinates of the (i+1)th waypoint are used to find the nearest grid point.In this way, AirTraf calculates the global fields of NO x and H 2 O emissions, fuel use and flight distance for output.After this step, the flight status check is performed at the end of the flying process.If the status is "arrived", the flight quits the flying process and its status is reset to "non-flight".Therefore, the flight status becomes either "inflight" or "non-flight" after the flying process.Once the status becomes "in-flight", the departure check is false in subsequent time steps t ≥ 2 and the aircraft moves to the new aircraft position without re-calculating the flight trajectory or fuel/emissions (Fig. 1, light blue).For simulations longer than 2 days, the same flight plan is reused: the departure time is automatically updated to the next day and the calculation procedures start from the departure check.

Fuel calculation
The methodologies of the fuel/emissions calculation module (Fig. 1, light orange) are described.Fuel use, NO x and H 2 O emissions are calculated along the flight trajectory obtained from the flight trajectory calculation.A total energy model based on the BADA methodology and the DLR fuel flow method is used.The fuel use calculation consists of the following two steps: a first rough trip fuel estimation and the second detailed fuel calculation (dashed boxes in Fig. 1, light orange).The former estimates an aircraft weight at the last waypoint (m n wp ), while the latter calculates fuel use for every flight segment and aircraft weights at any waypoint by backward calculation along the flight trajectory, using the m n wp as initial condition.
First, trip fuel (FUEL trip ) required for a flight between a given airport pair is roughly estimated: where FT is the estimated flight time (Table 2) and F BADA is the fuel flow.The BADA performance table provides cruise fuel flow data at specified flight altitudes for three different weights (low, nominal and high) under international standard atmosphere conditions.Hence, F BADA is calculated by interpolating the BADA data (assuming nominal weight) to the mean altitude of the flight (h, Table 2).Next, m n wp is estimated by where OEW, MPL and OLF are given in Table 1.The last term represents the sum of an alternate fuel, reserve fuel and extra fuel.It is assumed to be 3 % of the FUEL trip (r fuel = 0.03).The burn-off fuel required to fly from i = 1 to i = n wp and contingency fuel are assumed to be consumed during the flight and hence they are not included in m n wp .While the 3 % estimation is probably not far from reality for long-range flights, it is worth noting that typical reserve fuel quantities may amount to higher values, depending on the exact flight route.Airlines have their own fuel strategy and information about actual onboard fuel quantities is generally unavailable.Second, the burn-off fuel is calculated for every flight segment and the aircraft weights are estimated at all waypoints (the contingency fuel is disregarded in AirTraf (version 1.0)).With the BADA total energy model (Revision 3.9), the rate of work done by forces acting on the aircraft is equated to the rate of increase in potential and kinetic energy: where Thr and D are thrust and drag forces, respectively.m is the aircraft weight, g is the gravity acceleration, h is the flight altitude and dh/ dt is the rate-of-climb (or descent).For a cruise flight phase, both altitude and speed changes are negligible.Hence, dh/ dt = 0 as well as dV TAS /dt = 0 is assumed in AirTraf (version 1.0) and Eq. ( 3) becomes the typical cruise equilibrium equation: Thr i = D i at waypoint i.To calculate Thr i , the D i is calculated: where C L,i and C D,i are lift and drag coefficients, respectively.The performance parameters (S, C D0 and C D2 ) are given in Table 1, ρ i is the air density (Table 2) and V TAS,i is calculated at every waypoint (Table 2).The bank angle ϕ i is assumed to be zero.The thrust-specific fuel consumption (TSFC) η i and the fuel flow of the aircraft F cr,i are then calculated assuming a cruise flight: where C f1 , C f2 and C fcr are given in Table 1.The fuel use in the ith flight segment (FUEL i ) is calculated as where ETO i at the ith waypoint (in Julian date) is converted into seconds by multiplying with seconds per day (SPD, Table 1).The FUEL i incorporates the tail/head winds effect on V ground through ETO.The relation between the FUEL i and the aircraft weight (m i ) is obtained regarding the ith and (i + 1)th waypoints: Given m n wp by Eq. ( 2), the fuel use for the last flight segment FUEL n wp −1 and the aircraft weight at the last but one waypoint m n wp −1 can be calculated.This calculation is performed iteratively in reverse order from the last to first waypoints using Eqs.(3) to (10).Finally, the aircraft weight at the first waypoint m 1 is obtained.

Emission calculation
NO x and H 2 O emissions are calculated after the fuel calculations.NO x emission under the actual flight conditions is calculated by the DLR fuel flow method (Deidewig et al., 1996).It depends on the engine type, the power setting of the engine and atmospheric conditions.The calculation procedure follows four steps.First, the reference fuel flow of an engine under sea level conditions, f ref,i , is calculated from the actual fuel flow at altitude, f a,i (= F cr,i /(number of engines); see Eq. 8): www.
where δ total,i and θ total,i are correction factors.P total (in Pa) and T total (in K) are the total pressure and total temperature at the engine air intake, respectively, and P 0 and T 0 are the corresponding values at sea level (Table 1).P total and T total are calculated as where P a,i (in Pa) and T a,i (in K) are the static pressure and temperature under actual flight conditions at the altitude h i (Table 2).Here, h i is the altitude of the ith waypoint above the sea level (the geopotential altitude is used to calculate h i ).
The cruise Mach number M is given in Table 1.Second, the reference emission index under sea level conditions, EINO x,ref,i , is calculated using the ICAO engine emissions databank (ICAO, 2005) and the calculated reference fuel flow, f ref,i (Eq.11).Four data pairs of reference fuel flows f ref , and corresponding EINO x,ref , are tabulated in the ICAO databank for a specific engine under sea level conditions.Therefore, EINO x,ref,i values, corresponding to f ref,i , are calculated by a least squares interpolation (secondorder).
Third, the emission index under actual flight conditions, EINO x,a,i , is calculated from the EINO x,ref,i : H c,i = e (−19.0(qi −0.00634)) , (17) where δ total,i and θ total,i are defined by Eqs. ( 12) and ( 13), respectively.H c,i is the humidity correction factor (dimensionless number) and q i (in kg(H 2 O) (kg(air)) −1 ) is the specific humidity at h i (the unit ft is used here).Finally, NO x and H 2 O emissions under actual flight conditions are calculated for the ith flight segment using the calculated FUEL i (Eq.9): where the Penner et al., 1999).The H 2 O emission is proportional to the fuel use, assuming an ideal combustion of jet fuel.The NO x and H 2 O emissions are included in the flight properties (Table 2).With regard to the reliability of the fuel/emissions calculation using these methods, Schulte et al. (1997) showed a comparison of measured and calculated EINO x for some aircraft/engine combinations (Schulte et al., 1997).The study gave some confidence in the prediction abilities of the DLR method, although it showed that the calculated values from the DLR method underestimated the measured values on average by 12 %.In Sect. 5 we verify the methods, using 1-day AirTraf simulation results.Detailed descriptions of the total energy model and the DLR fuel flow method can be found elsewhere (Eurocontrol, 2011;Deidewig et al., 1996).

Aircraft routing methodologies
The current aircraft routing module (Fig. 1, light green) works only with respect to the great circle and flight time routing options.These routing methodologies are described in Sects.3.1 and 3.2.Benchmark tests are performed offline (without EMAC) to verify the accuracy of the methodologies.

Formulation of great circles
AirTraf calculates a great circle at any arbitrary flight altitude with the great circle routing option.First, the coordinates of the waypoints are calculated.For the ith and (i + 1)th waypoints, the central angle σi (Vincenty, 1975): where φ i (in rad) is the latitude of the ith waypoint and λ i (in rad) is the difference in longitude between the ith and (i + 1)th waypoints.The Vincenty formula was set as the default method, while optionally the spherical law of cosines or the Haversine formula can be used in AirTraf to calculate σ (unshown).With Eq. ( 21), the great circle distance for the ith flight segment d i is calculated: or For the great circle routing option, flight altitudes at all waypoints are set as h i = constant for i = 1, 2, • • •, n wp (h i is used in kilometers in Eqs.22 and 23) and either Eq. ( 22) or ( 23) is used to calculate d i .Equation ( 22) calculates d i by an arc and hence the great circle distance between airports, i.e., n wp −1 i=1 d i , is independent of n wp .On the other hand, Eq. ( 23) calculates d i by linear interpolation in polar coordinates.In that case, n wp −1 i=1 d i depends on n wp ; the sum becomes close to that calculated from Eq. ( 22 3. are compared to those with other routing options, Eq. ( 23) should be used for the comparison with the same n wp .In addition, Eq. ( 23) is used for the flight trajectory optimization (see Sect. 3.2), because it is necessary to calculate d i including altitude changes.
Next, the true air speed V TAS and the ground speed V ground at the ith waypoint are calculated: where M is the Mach number, γ is the adiabatic gas constant and R is the gas constant for dry air (Table 1).Temperature T i and three-dimensional wind components (u i , v i , w i ) of the ith waypoint are available from the EMAC model fields at t = 1; the local speed of sound a i is then calculated (Table 2).The flight direction is calculated for every flight segment by using the three-dimensional coordinates of the ith and (i + 1)th waypoints.Thereafter, V TAS,i , V wind,i and V ground,i (scalar values) corresponding to the flight direction are calculated.As shown in Eq. ( 25), the influence of tail/head winds on ground speed is considered.In AirTraf, M was set constant as default.It is also possible to perform AirTraf simulations with different options, such as V TAS,i = constant and V wind,i = 0. Finally, ETO i (in Julian date) and FT (in s) are calculated as where ETO 1 is the departure time of the flight and ETO i incorporates the influence of tail/head winds on the flight.

Benchmark test on great circle calculations
A benchmark test of the great circle routing option was performed to confirm the accuracy of the great circle distance calculation.Great circles were calculated for five representative routes without EMAC (offline).Table 3 shows the information for the five routes (the locations are shown in Fig. 4).The characteristics of the routes were as follows: R1 consisted of an airport pair in the Northern Hemisphere (MUC-JFK) and the difference in longitude between them was λ airport < 180 • ; R2 consisted of an airport pair in the Northern Hemisphere (HND-JFK) with λ airport > 180 • (discontinuous longitude values due to the definition of the longitude range [−180, 180]); R3 consisted of an airport pair in the Northern and Southern Hemisphere (MUC-SYD); R4 was a special route, where λ airport = 0 • and the difference in latitude was φ airport = 0 • ; and R5 was another special route with λ airport = 0 • and φ airport = 0 • .Other calculation conditions were set as follows: M = 0.82, h i = 0, a i = 304.5 ms −1 and V TAS,i = V ground,i = 249.7 ms −1 (under nowind conditions, i.e., V wind,i = 0) for i = 1, 2,  22) and ( 23), and were compared to that calculated with MTS.In addition, the sensitivity of the great circle distance with respect to n wp was analyzed varying n wp in the range [2, 100].
Table 4 shows the calculated great circle distances by Eqs. ( 22) and ( 23) and MTS.The columns 5 to 7 show the difference in the distance among them (see the caption of Table 4 for more details).The results showed that both d eq23,eq22 and d eq23,MTS varied between −0.0036 and −0.0005 %, while d eq22,MTS showed 0.0 %.The great circle distances calculated by Eqs. ( 22) and ( 23) were accurate to −0.004 %, and hence this routing option works properly.
Figure 5 shows the result of the sensitivity analysis of n wp on the great circle distance.The results show that the distance calculated by Eq. ( 22) (open circle) has no dependence on n wp as noted in Sect.3.1.1,whereas that by Eq. ( 23) (closed circle) depends on n wp and converged with increasing n wp : the accuracy of the results by Eq. ( 23) decreased when using fewer n wp .For n wp ≥ 20, the results of Eqs. ( 22) and ( 23) were almost the same.Therefore, n wp ≥ 20 is practically desired for the use of Eq. ( 23).

Overview of the genetic algorithm
The flight trajectory optimization with respect to the flight time was performed using GA (Holand, 1975;Goldberg, 1989), which is a stochastic optimization algorithm.The Aircraft routing module (Fig. 1, light green) is linked to the flight trajectory optimization module (Fig. 1, dark green); this optimization module consists of the Adaptive Range Multi-Objective Genetic Algorithm (ARMOGA version 1.2.0) developed by D. Sasaki and S. Obayashi (Sasaki et al., 2002;Sasaki andObayashi, 2004, 2005).ARMOGA will be implemented as part of the MESSy infrastructure in the next version of MESSy so that it can be used for optimization problems by other submodels as well.For each routing option (except for the great circle routing option), a single-objective optimization problem is solved.The main advantage of GA shows the result calculated with the Movable Type scripts (MTS), using the Haversine formula with a spherical Earth radius of R E = 6371 km.
In columns 5 to 7:  is that GA requires neither the computation of derivatives or gradients of functions, nor the continuity of functions.Therefore, various evaluation functions (called objective functions) can easily be adapted to GA.As for the working principle of GA, a random initial population is created and the population evolves over generations to adapt to an environment by the genetic operators: evaluation, selection, crossover and mutation.When this biological evolutionary concept is applied for design optimizations, fitness, individuals and genes correspond to an objective function, solutions and design vari-ables, respectively.A solution found in GA is called an optimal solution, whereas a solution with the theoretical optimum of the objective function is called the true-optimal solution.If GA works properly, it is expected that the optimal solution will converge to the true-optimal solution.On the other hand, the main disadvantage of GA is that GA is computationally expensive.The flight trajectory optimization is applied for all flights and therefore a user has to choose appropriate GA parameter settings to reduce computational costs (or find a compromise for the settings, which sometimes depend on the computing environment).

Formulation of flight trajectory optimization
The flight trajectory optimization is described focusing on geometry definitions of the flight trajectory, the definition of the objective function and the genetic operators.There exist a number of selection, crossover and mutation operators in ARMOGA.Therefore, the genetic operators employed in this study are described here.
A solution x (the term is used interchangeably with flight trajectory) is a vector of n dv design variables: x = (x 1 , x 2 , • • •, x n dv ) T .Using the design variable index j (j = 1, 2, • • •, n dv ), the j th design variable varies in lower/upper bounds [x l j , x u j ].GA searches for the optimal solution, corresponding to the routing option, around the great circle of an airport pair including altitude changes.Figure 6 shows the geometry definition of a flight trajectory from MUC to JFK as an example: the projection on the Earth (bottom) with three control points (CPs, black circles) and the vertical cross section (top) with five CPs.The coordinates of the airports were given from a flight plan (Fig. 1, dark blue) and were  6.Bottom: the dashed boxes show rectangular domains of three control points.: central points of the domains are calculated on the great circle (thin solid line), which divide the λ airport into four equal parts.Top: the dashed lines show the lower/upper variable bounds in altitude."FL290" stands for a flight level at 29 000 ft. Longitude coordinates for x 7 , x 8 , • • •, x 11 are pre-calculated; the coordinates divide the λ airport into six equal parts.fixed (the coordinates of MUC and JFK are shown in Table 5).
Six design variables x j (j = 1, 2, • • •, 6) were used for location, as shown in Fig. 6 (bottom).x 1 , x 3 and x 5 indicate longitudes, while x 2 , x 4 and x 6 indicate latitudes.To create three rectangular domains for the design variables (dashed boxes), central points of the domains (diamond symbols) were calculated.The points are located on the great circle, dividing the longitude distance between MUC and JFK ( λ airport ) into four equal parts.After that, the three domains centered around the central points were created.The domain size was set to 0.1 × λ airport (short-side) and 0.3 × λ airport (longside).This procedure calculates the lower/upper bounds of the six design variables, i.e., [x l j , x u j ] (j = 1, 2, • • •, 6), and Table 6 lists these values.GA provided the values for x 1 to x 6 within the respective bounds (i.e., the values were generated within the rectangular domains) and the coordinates of the three CPs were determined: CP1 (x 1 , x 2 ), CP2 (x 3 , x 4 ) and CP3 (x 5 , x 6 ).A flight trajectory is represented by a B-spline curve (third-order) with the three CPs as locations (bold solid line, Fig. 6 bottom), and then any arbitrary number of waypoints is generated along the trajectory.To generate the same number of waypoints between the CPs, n wp was calculated as mod(n wp − 1, n CP loc + 1) = 0, where the number of CPs was n CP loc = 3.
For the altitude direction, five design variables x j (j = 7, 8, • • •, 11) were used (Fig. 6, top).Here x 7 to x 11 indicate altitude values.With the lower h l and the upper h u variable bound parameters, the bounds of the five design variables were determined by x l j = h l and x u j = h u for j = 7, 8, • • •, 11.In this study, h l = FL290 and h u = FL410, as listed in Table 6 ("FL290" stands for a flight level at 29 000 ft).These altitudes correspond to a general cruise flight altitude range of commercial aircraft (Sridhar et al., 2013).GA provided the values of x 7 to x 11 in [FL290, FL410] and the coordinates of the five CPs were determined: CP4 (x 7 ), CP5 (x 8 ), CP6 (x 9 ), CP7 (x 10 ) and CP8 (x 11 ).Note that these values vary freely between FL290 and FL410 to explore widely the possibility of minimizing climate impact by aircraft routing.The longitude coordinates of the five CPs were pre-calculated to divide the λ airport into six equal parts.The altitudes of the airports were fixed at h l (= FL290).A flight trajectory is also represented by a B-spline curve (third-order) with the five CPs in the vertical cross section (bold solid line, Fig. 6 top) and then waypoints are generated along the trajectory in such a way that the longitude of the waypoints is the same as that for the flight trajectory projected on the Earth.
GA starts its search with a random set of solutions (population approach).The initial population operator (Fig. 1, dark green) provides initial values of the 11 design variables at random within the lower/upper bounds described above, thereby creating solutions.The operator creates n p different solutions (where n p is the population size).To evaluate the solutions, the objective function f was calculated for each of the solutions by summing the flight time over all flight segments (Fig. 1, dark green).The single-objective optimization problem on the flight time can be written as follows: where n dv = 11, d i and V ground,i are calculated by Eqs. ( 23) and ( 25), respectively (V TAS,i and V wind,i are calculated as described in Sect.3.1.1).No constraint function is used in AirTraf (version 1.0).Good solutions are identified in the population by Fonseca and Fleming's Pareto ranking method (Fonseca and Fleming, 1993), although the single-objective optimization is solved here.A rank of a solution was assigned proportionally to the number of solutions that dominate it, and a fitness value of a solution was computed by 1/ rank (no fitness sharing was used).A solution with a higher fitness value (i.e., a smaller rank value) has a higher probability of being copied into a mating pool.The stochastic universal sampling selection (Baker, 1985) makes duplicates of good solutions in the mating pool at the expense of bad solutions based on cumulative probability values, while keeping the size of n p .
x j,c1 and x j,c2 denote the j th design variable of the child solutions, and x j,p1 and x j,p2 denote the j th design variables of the parent solutions (the mated pair of the old generation).α is a user-specified crossover parameter and u 1 is a random number between zero and one.Thereafter, the mutation operator added a disturbance to the child solutions by the revised polynomial mutation operator (Deb and Agrawal, 1999) with a mutation rate r m .A polynomial probability distribution was used and the mutated design variable was created.The parameter δ q is first calculated as where δ = min[(x j,c −x l j ), (x u j −x j,c )]/(x u j −x l j ).The j th design variable varies in [x l j , x u j ]. u 2 is a random number be-tween zero and one, and η m is an external parameter controlling the shape of the probability distribution.The mutated design variable (mutated child solution) x j,mc is calculated as follows: Using the genetic operators above, it is expected that the population of solutions will be improved and a new and better population created in subsequent generations.When the evolution is computed for a fixed number of generations n g , GA quits the optimization and an optimal solution showing the best f of the whole generation is output.The optimal solution has the superior combination of the 11 design variables x = (x 1 , x 2 , • • •, x 11 ) T to minimize f .The flight properties of the optimal solution are also available (ETO, h, FT, etc., listed in the first and second groups (divided by rows) of Table 2).The flight trajectory optimization methodology described here could be applied to any routing option (except for the great circle routing option).In that case, the objective function f given by Eq. ( 28) needs to be reformulated corresponding to the selected routing option.

Benchmark test on flight trajectory optimization with flight time routing option
To quantify the performance of GA, there is a need to choose an appropriate benchmark test of the flight trajectory optimization, where the true-optimal solution f true of the test is known.Here, the single-objective optimization for minimization of flight time from MUC to JFK was solved without EMAC (offline); that is, the optimization problem defined in Sect.3.2.2 was solved.Calculation conditions for the test are summarized in Table 5. V wind was set to 0 km h −1 (nowind conditions); V TAS and V ground were set to 898.8 km h −1 (constant).Hence, f true equals the flight time along the great circle from MUC to JFK at FL290 (having its minimum d i in

Optimization results
The influence of the population size n p and the number of generations n g on the convergence properties of GA was examined.Figure 7 shows the optimal solutions varying with n g for a number of fixed n p .The results confirmed that the optimal solutions come sufficiently close to f true with increasing n p and n g .The optimal solution showing the closest flight time to f true was obtained for n p = 100 and n g = 100.This solution is called best solution in this study and its flight time was f best = 25 996.6 s.The difference in flight time between the f best and f true was f < 3.0 s (less than 0.01 %).
To confirm the diversity of GA optimization, we focus on the optimization yielding the best solution (n p = 100 and n g = 100).Figure 8 shows all the solutions explored by GA.It is clear that GA explored diverse solutions from MUC to JFK including altitude changes and found the best solution.As shown in Fig. 8, the best solution (red line) overlapped with the true-optimal solution, i.e., the great circle at FL290 (dashed line, black).To investigate the difference between the solutions, the comparisons of trajectories for the best solution and the true-optimal solution in the vertical cross section are plotted in Fig. 9.The maximum difference in altitude is less than 1 m.Therefore, GA is adequate for finding an op- timal solution with sufficient accuracy (in a strict sense, this conclusion is confined to the benchmark test).

Dependence of initial populations
To analyze the dependence of the optimal solution on the initial population, Fig. 10    of objective function evaluations (= n p × n g ) for the 10 independent GA simulations from different initial populations with n p = 100 and n g = 100.Figure 10 shows that the 10 solutions converged in early generations and gradually continued to converge to f true with an increasing number of function evaluations.The convergence behavior is similar among the 10 simulations, regardless of the initial population.Ta-  S1 summarizes the 10 optimal solutions in detail.
ble S1 in the Supplement shows a summary of the 10 optimal solutions.As indicated in Table S1, the value of the objective function f (flight time) is slightly different.f (= f − f true ) ranged from 2.5 to 3.7 s, which is approximately 0.01 % of f true .In addition, the mean value of the 10 objective functions was f = 2.9 s (0.01 % of f true ) and the standard deviation was s f = 0.4 s (0.001 % of f true ).Therefore, the variation in the objective function with different initial populations is small.

Population and generation sizing
With increased n p and n g , GA tends to find an improved solution.It is important to note that the required size of n p and n g is problem-dependent.However, following a simple initial guess for n p and n g is a good starting point for their sizing.The influences of n p and n g on the accuracy of GA optimizations and on the variation in the optimal solution due to different initial populations were analyzed.Figure 11 shows the f and s f for all the combinations of n p and n g .The results confirm that f and s f decrease with an increase in n p and n g .That is, the optimal solution converges to the trueoptimal solution (the accuracy increases) and the variation in the optimal solution due to different initial populations decreases (the dependency decreases).
On the other hand, computational costs should also be kept as low as possible for practical use of EMAC/AirTraf (online) applied to long-term global air traffic simulations.Figure 12 shows the variation of f and s f for all combina- , where n = 10.f and s f (in %) relative to the true-optimal solution are calculated as ( f /f true )×100 and (s f /f true )×100, respectively.tions of n p and n g with respect to the number of function evaluations.The symbols and error bars in the figure correspond to f and s f , respectively (Table S2 in the Supplement lists these values).The results showed that there is a trade-off between the accuracy of GA optimizations and the number of function evaluations (i.e., computing time).The figure also shows the power function (red line) fitted to the results by using the standard least squares algorithm (see the caption in Fig. 12 for more details).The enlarged drawing in Fig. 12 shows that if one selects the number of function evaluations (= n p × n g ) of 800, the large reduction of computational costs of 92 % can be achieved, keeping f less than 0.05 % (s f ≈ 0.02 %), compared to the optimal solution obtained by 10 000 function evaluations (n p = 100 and n g = 100).For n p × n g = 800, one can select any combination of n p and n g : n p = 10 and n g = 80, n p = 20 and n g = 40, etc.A user makes his/her own choice on n p and n g by referring to the values of f and s f shown in Fig. 12.Similarly, a reduction of 97 % can be achieved, keeping f less than 0.1 % (s f ≈ 0.04 %).Therefore, computational costs can be reduced drastically by selecting n p and n g for different purposes.

Demonstration of a 1-day AirTraf simulation
The aircraft routing methodologies corresponding to the great circle and flight time routing options were verified in Sect.3. Here, 1-day AirTraf simulations were performed in EMAC (online) with the respective routing options for demonstration.
Figure 12.Chart for finding the appropriate number of function evaluations (= n p ×n g ), including the enlarged drawing in the early 1500 evaluations.The symbols with error bars correspond to f ± s f (in %); their definitions are given in the caption in Fig. 11.The fitted curve (power function, red line) to f is y = e 0.92 x −0.59 , where x are the function evaluations and y is f (in %); R 2 = 0.89.The fitted curve to s f is calculated similarly: y = e 0.67 x −0.73 , where R 2 = 0.71 (unshown).

Simulation setup
We focus on the trans-Atlantic region for the demonstration, because the optimization potential is possibly large for this region.Table 7 lists the setup for the 1-day simulations.The simulations were performed for 1 typical winter day in the T42L31ECMWF resolution.The weather situation on that day showed a typical weather pattern for winter characterized by westerly jet streams in the North Atlantic region.
The number of trans-Atlantic flights in the region was 103 (52 eastbound flights and 51 westbound flights).We assumed that all flights were operated by A330-301 aircraft with CF6-80E1A2 (2GE051) engines.Thus, the data shown in Table 1 were used.Four 1-day simulations were separately performed for the great circle routing option at fixed altitudes FL290, FL330, FL370 and FL410 (see Sect. 3.1.1).In addition, a single 1-day simulation was performed for the flight time routing option, including altitude changes in the range of [FL290, FL410] (see Sect. 3.2.2).For the two options, the Mach number was set to M = 0.82 and therefore the values of V TAS and V ground were different at every waypoint (Eqs.24 and 25).The number of waypoints was set to n wp = 101.As described in Sect.3.1.1,the flight distance was calculated by Eq. ( 23) for the two routing options.The optimization parameters were set as follows: n p = 100, n g = 100 and other GA parameters were the same as those used in the benchmark test in Sect.3.2.3.
www.geosci-model-dev.net/9/3363/2016/Geosci.Model Dev., 9, 3363-3392, 2016 Table 8.Information for the trajectories of the three selected airport pairs; they were extracted from the 1-day AirTraf simulations.Columns 7 to 11 show the obtained flight times for the flight time and great circle routing options.GC FL290: great circle at 29 000 ft.The 1-day simulation was parallelized on four PEs of Fujitsu Esprimo P900 (Intel Core i5-2500 CPU with 3.30 GHz; 4 GB of memory; peak performance of 105.6 × 4 GFLOPS) at the Institute of Atmospheric Physics, German Aerospace Center.The 1-day simulation required approximately 15 min for the great circle routing option, while it took approximately 20 h for the flight time routing option.Most of the computational time is consumed by the trajectory optimizations.Therefore this time can be reduced by choosing properly all GA parameters, using more PEs, or decreasing n p and n g .As discussed in Sect.3.2.6, a large reduction in computing time of roughly 90 % can be achieved by using a small n p and n g with still sufficient accuracy of the optimizations.

Optimal solutions for selected airport pairs
The 1-day simulation results for the flight time routing option confirmed that the optimized flight trajectories showed a large altitude variation.To give an overview of the optimizations, we classified those optimized flight trajectories according to their altitude changes into three categories.Type I: eastbound and westbound time-optimal flight trajectories showed little altitude changes; Type II: the eastbound time-optimal flight trajectory showed little altitude changes, while the westbound time-optimal flight trajectory showed distinct altitude changes; and Type III: eastbound and westbound time-optimal flight trajectories showed distinct altitude changes.We have selected three airport pairs of each type and Table 8 shows the details of them.Here, we mainly discuss the selected solutions of Type II, which were eastbound and westbound flights between Minneapolis (MSP) and Amsterdam (AMS).
We examined first the optimal flight trajectories between MSP and AMS. Figure 13 shows all trajectories explored by GA (black lines) and the time-optimal flight trajectories for eastbound and westbound flights (red and blue lines).Figure 13a and b show that GA explored diverse trajectories properly considering altitude changes in the range of [FL290, FL410].Similar results were obtained for the selected solutions of Types I and III, as shown in Figs.S1 and S2 in the Supplement.In addition, the eastbound time-optimal flight trajectory was located at FL290, while that for westbound flights showed large altitude changes; i.e., it climbed, descended, climbed and then descended again.The mean flight altitudes of these trajectories were h = 8839 m and h = 10 002 m.These time-optimal flight trajectories were compared to the prevailing wind fields.To calculate tail/head winds in the eastern and western directions, the major wind component is shown in Fig. 14.The contours represent the zonal wind speed (u); black arrows show the wind speed ( √ u 2 + v 2 ) and direction at the departure time at h. Figure 14a and b show that the eastbound time-optimal flight trajectory (red line) was located to the south of the great circle (black line) to take advantage from the tail winds of the westerly jet stream (red region), while the westbound time- optimal flight trajectory (blue line) was located to the north of the great circle to avoid the head winds (red region).Similar comparisons for the selected solutions of Types I and III showed that the obtained optimal flight trajectories effectively take advantage of the wind fields (see Supplement, Figs.S3 and S4).
To understand the behavior of the altitude changes of the optimal flight trajectories, Fig. 15 shows the altitude distribution of the true air speed (V TAS ) and the tail wind indicator (V ground /V TAS ) along the time-optimal flight trajectories.The indicator was calculated by Eq. ( 25) transformed into V ground /V TAS = 1 + V wind /V TAS ; this means tail winds ((V ground /V TAS ) ≥ 1.0) and head winds ((V ground /V TAS ) <  1.0) in the flight direction.Figure 15c shows that the core tail winds region was located at 8.5 km and the tail winds were most beneficial for the eastbound flight trajectory.On the other hand, the westbound flight trajectory went through the regions where V TAS was high, as shown in Fig. 15b.In addition, Fig. 15d shows that the descent at a flight time of 16 000 s was effective to counteract the head winds.These results confirm that GA correctly takes into account the weather conditions and finds the appropriate flight trajectories corresponding to the flight direction.Similar results were www.geosci-model-dev.net/9/3363/2016/Geosci.Model Dev., 9, 3363-3392, 2016 Figure 16.Flight time (in %) vs. number of function evaluations (= n p × n g ) for three selected airport pairs, including the enlarged drawing in the early 1500 evaluations.The population size n p is 100 and the number of generations n g is 100.f * means the difference in flight time between the solution f and the obtained optimal solution f opt , which was finally obtained after 10 000 function evaluations.This was chosen because f true for the six flights are unknown.The f opt for each flight corresponds to the flight time for the timeoptimal case (column 7, Table 8).The f * (in %) is calculated as obtained for the solutions of Types I and III (see Supplement, Figs.S5 and S6).
Next, we compared the resulting flight times for the selected solutions.Table 8 shows the obtained flight times for the time-optimal and great circle cases.As shown in Table 8, the flight time is lower for the time-optimal case compared to the great circle cases.In addition, the flight time is lower for the eastbound time-optimal flight trajectories compared to that for the westbound time-optimal flight trajectories.This supports the observation that GA correctly takes into account weather conditions for the trajectory optimization.With regard to the convergence behavior of the optimization, Fig. 16 shows the flight time vs. the number of objective function evaluations corresponding to the GA simulations for the three selected airport pairs.As expected, the solutions converged to each optimal solution.Thus, GA successfully found the time-optimal flight trajectories for the three airport pairs.It is also clear from Fig. 16 that a reduction in computing time can be achieved by choosing properly n p and n g , although the solutions converged more slowly under the wind conditions than those under no-wind conditions (Fig. 12).

One-day simulation results for all flights
Next, the 1-day simulation results for 103 trans-Atlantic flights are analyzed.Figure 17 shows the obtained flight tra-jectories for the flight time and great circle routing options.Figure 17a and c show that many eastbound time-optimal flight trajectories congregated around 50 • N over the Atlantic Ocean to take advantage from the tail winds in the westerly jet stream.On the other hand, the westbound time-optimal flight trajectories were located to the north and south of that region to avoid head winds (as shown in Fig. 17b and d).
In addition, Fig. 17a and b show that only 5 of 52 eastbound time-optimal flight trajectories showed large altitude changes, in comparison to 35 of 51 westbound time-optimal flight trajectories.The mean flight altitudes for the 52 eastbound, 51 westbound and total 103 flights were h = 9029 m, 9517 m and 9271 m, respectively.
As shown in Fig. 15, altitude changes were due to variations of V TAS and prevailing winds.We now confirm this behavior, focusing on the results for all flights.Figure 18a and  b show the values of V TAS and V ground /V TAS at waypoints for the time-optimal and great circle flights, with linear lines fitted by the least squares algorithm.Figure 18a shows that V TAS is higher at low altitudes.From Eq. ( 25), high V TAS values increase V ground values, thereby minimizing flight time.The mean V TAS for the time-optimal and great circle cases are shown in Table 9.The mean V TAS value (column 4) for the time-optimal case is 245.1 ms −1 , while that for the great circle cases ranges from 241.2 to 244.9 ms −1 , although the mean flight altitude for the time-optimal case is h = 9271 m, which is higher than FL290 (= 8839 m).GA successfully found the flight trajectories with high V TAS values as timeoptimal flights.
With regard to the wind effects, Fig. 18b shows that the fitted line for the eastbound time-optimal case (solid line, red) is larger between FL290 (= 8839 m) and 9500 m compared to that for the eastbound great circle case (dashed line, red).These altitude bounds are effective under the present weather condition to take advantage of tail winds for the eastbound flights.Thus, almost all the eastbound time-optimal flight trajectories were located at FL290, as shown in Fig. 17a (top).On the other hand, the fitted line for the westbound timeoptimal case (solid line, blue) is distributed widely in altitude and is larger between FL290 (= 8839 m) and 12 000 m compared to that for the westbound great circle case (dashed line, blue).The westbound time-optimal flight trajectories certainly mitigated the head winds effect.Thus, many westbound time-optimal flight trajectories showed large altitude changes, as shown in Fig. 17b (top).The similar plot of V ground is shown in the Supplement (Fig. S7), which incorporates the influences of both V TAS and winds; the plot indicates similar trends as shown in Fig. 18b.Table 9 also shows that the mean V ground value (column 7) for the timeoptimal case is 250.2 ms −1 , while that for the great circle cases ranges from 241.1 to 244.7 ms −1 .Therefore, the trajectories found by GA through altitude changes passed areas, which correctly lead to larger V ground .
These altitude changes affect the fuel consumption (the term is used interchangeably with fuel flow).Figure 19 shows Table 9.The mean value of V TAS and V ground for the time-optimal and great circle cases.The mean values were calculated using V TAS and V ground values at all waypoints.Eastbound: average of 52 eastbound flights; Westbound: average of 51 westbound flights; and Total: average of 103 flights.
Table 10.The mean fuel consumption (in kg(fuel)min −1 ) for the time-optimal and great circle cases.Eastbound: average of 52 eastbound flights; Westbound: average of 51 westbound flights; and Total: average of 103 flights.Columns 5 to 7 show the reference cruise fuel consumption (in kg(fuel)min −1 ) for three different weights (low, nominal and high) in the international standard atmosphere.BADA provides the reference data at specific flight altitudes.Therefore, the reference values for the time-optimal case in parentheses were estimated from the reference data at FL290 and FL330 by linear interpolation (the mean flight altitude of the time-optimal case was h = 9271 m, which is the value between FL290 (= 8839 m) and FL330 (= 10 058 m)).Linear fits of the time-optimal (solid line, red (eastbound) and blue (westbound)) and great circle cases (dashed line, red (eastbound) and blue (westbound)) are included.V TAS of the international standard atmosphere (ISA) is given in (a) (solid line, black) provided by the BADA atmosphere table (Eurocontrol, 2010).

Case
sults show that all symbols lay on the right side of the 1 : 1 solid line.That is, the flight time for the time-optimal flights is lower compared to that for the great circle flights for all airport pairs.Table 11 shows the total flight time simulated by AirTraf for eastbound, westbound and total flights.The total value is certainly minimal for the time-optimal case, while in relative terms the value increases by +1.5, +2.5, +2.9 and +2.9 % for the great circle cases at FL290, FL330, FL370 and FL410, respectively.Regarding the total value of fuel use, Table 11 indicates that the value increases by +5.4 % for the great circle case at FL290 when compared with the value of the time-optimal case.On the other hand, the fuel use decreased by −5.8, −14.9 and −20.8 % for the great circle cases at FL330, FL370 and FL410, respectively.The total values of NO x and H 2 O emissions show a similar trend:  10. the total value of NO x emission increased by +5.2 % for the great circle at FL290, while it decreased by −12.9, −24.9 and −29.4 % for the great circle cases at FL330, FL370 and FL410, respectively.The changes in total H 2 O emission were the same as those of the total fuel use, because EIH 2 O = 1230 g(H 2 O) (kg(fuel)) −1 was used.Figure 19 already shows that the mean fuel consumption for the timeoptimal case is high, owing to the low mean flight altitude.Thus, the total amount of fuel use increased for this case, which increased total NO x and H 2 O emissions.and H 2 O emissions are not representative for all seasons and the whole world's air traffic, because they have been obtained under the specific winter conditions using the trans-Atlantic flight plans.

Verification of the AirTraf simulations
To verify the consistency for AirTraf simulations, the 1-day simulation results described in Sect. 4 are compared to reference data of flight time, fuel consumption, EINO x and aircraft weight.Data obtained under similar conditions (aircraft/engine types, flight conditions, weather situations, etc.) were selected for the comparison, although the conditions are not completely the same as the calculation conditions for the 1-day simulations.Note that the verification of the aircraft weight is related to that of the fuel use calculations, because the aircraft weight was calculated by adding the amount of fuel use (Eq.10).In addition, H 2 O emission is proportional to the fuel use assuming ideal combustion.Thus, its verification would be redundant.First, Table 12 shows the flight time for the seven timeoptimal flight trajectories simulated by AirTraf and three reference data (the seven airport pairs are geographically close to those of the reference data).Sridhar et al. (2014) simulated the wind-optimal flight trajectory from Newark (EWR) to Frankfurt (FRA) using a specific winter day, and the flight time was 22 980 s.The flight time of the time-optimal flight trajectory from JFK to FRA simulated by AirTraf was 22 955 s.This agrees well with the value reported by Sridhar et al. (2014).Irvine et al. (2013) analyzed the variation in flight time of time-optimal flight trajectories between JFK and London (LHR) using weather data for three winters.The results showed that the flight time for eastbound and westbound flights ranged from approximately 18 000 to 22 200 s, and from 21 600 to 27 000 s, respectively (see Fig. 3 in Irvine et al., 2013).In addition, Grewe et al. (2014a) optimized the trans-Atlantic 1-day air traffic (for winter) with respect to air traffic climate impacts and economic costs to investigate routing options for minimizing the impacts.The results showed that the mean flight time of the air traffic ranged from 26 136 to 27 792 s (eastbound), while it ranged from 29 664 to 31 788 s (westbound), depending on the degree of climate impact reduction (see Tables 2 and 3 in Grewe et al., 2014a).The flight times between the seven airport pairs are close to the reference data and the variation shows a good agreement with the trend of the increased flight times for westbound trans-Atlantic flights in winter due to westerly jet streams, as indicated from the reference data.
Second, the fuel consumption was verified using the mean fuel consumption value of 103 flights and the reference data, as shown in columns 4 to 7 of Table 10.Note that the AirTraf simulations were performed under the specific winter conditions (Table 7), while the reference data show the estimated values under international standard atmosphere conditions.Table 12.The flight time for time-optimal flight trajectories from 1-day AirTraf simulations and optimal trajectories from earlier studies.The original units of the flight times of the studies are converted into seconds.
Table 10 shows that the mean fuel consumption values for the time-optimal and great circle cases (column 4) were comparable to those of the reference data corresponding to low and nominal weights (columns 5 and 6).In the AirTraf simulations, the overall load factor of the worldwide air traffic was used (Table 1).If a specific load factor of A330-301 for international flights is available, the value is possibly higher than 0.62 and the corresponding mean fuel consumption values are expected to increase.
Third, the mean EINO x (in g(NO x ) (kg(fuel)) −1 ) simulated by AirTraf were compared to the six reference data.Table 13 shows that the obtained mean EINO x value is lower at high altitudes, and it ranged from 10.8 to 12.2 g(NO x ) (kg(fuel)) −1 .These values are in the same range as the reference data.Note that the reference data provided by Sutkus et al. (2001) show higher EINO x values.They correspond to the values for the CF6-80E1A2 (1GE033) engine instead of the CF6-80E1A2 (2GE051) engine used in our simulations.NO x of aircraft engines, in general, decreases owing to an installation of a new combustor.The 2GE051 utilizes the new 1862M39 combustor, which is known as a low-emissions combustor.Thus, the reference EINO x value of 2GE051 will be lower than that of the 1GE033.
Finally, the aircraft weights simulated by AirTraf were verified to make sure that the fuel use calculations were performed properly.AirTraf simulates realistic fuel consumptions under cruising flight; i.e., the aircraft weight decreases from the first waypoint (m 1 ) to the last waypoint (m n wp ) as fuel is burnt (as described in Sect.2.5).Thus, m 1 and m n wp correspond to the maximum and minimum aircraft weights, respectively.Here the obtained m 1 and m n wp for the 103 flights were compared with three structural weight limits (MTOW, MLW and MZFW), which are commonly used to provide flight operations safety, and one specified weight limit (MLOW) of the A330-301 aircraft.Table 14 shows the designated constraints among the m 1 , m n wp and the four weight limits.Note that no model that constrains the structural weight limits was included in AirTraf.
As indicated in Table 14, the first constraint is on maximum take-off weight (MTOW).The MTOW is limited for the aircraft so as not to cause structural damage to the airframe during take-off.Figure 21 shows a comparison of m 1 and m n wp with the weight limits (MTOW, MLW and MLOW).The results showed that almost all the m 1 (closed circles) were less than the MTOW.Only 15 of 515 flights (total of the time-optimal and great circle cases: 5 cases × 103 flights) exceeded the MTOW.For these 15 flights, actual flight planning data indicate higher flight altitudes to increase the fuel mileage, leading to the decrease in m 1 .The second constraint is on maximum landing weight (MLW).To prevent structural damage to the landing gear and the fuselage, an aircraft has to reduce the total weight until MLW prior to landing. Figure 21 shows that all the m n wp (open circles) were certainly less than MLW.The third constraint is on maximum zero fuel weight (MZFW), which corresponds to  14.
the maximum operational weight of the aircraft without usable fuel.The MZFW of an A330-301 aircraft is 164 000 kg (EASA, 2013), while the calculated zero fuel weight (ZFW) was 154 798 kg for all flights.This always satisfies the third constraint ZFW ≤ MZFW.Note that the ZFW is calculated as ZFW = OEW + MPL × OLF, and hence it depends only on the aircraft type and the load factor (Table 1).In addition, the fourth constraint is on the approximately minimum operational weight of an A330-301 aircraft in the international standard atmosphere (MLOW).The MLOW is used here as a measure of validity of fuel use calculations and is not a strict constraint.As shown in Fig. 21, all the m n wp (open circles) were higher than the MLOW.As a result, almost all the m 1 and m n wp simulated by AirTraf satisfied the four constraints.Thus, AirTraf simulates fairly good fuel use.

Conclusions
This study presents the AirTraf (version 1.0) global air traffic submodel of EMAC.The great circle and flight time routing options can be used in AirTraf 1.0.Two benchmark tests were performed without EMAC (offline).First, a benchmark test was performed for the great circle routing option using five representative routes.The results showed that the routing methodology works properly and the great circle distances showed quantitatively good agreement with those calculated by MTS.The accuracy of the results was within −0.004 %.Second, a benchmark test was performed for the flight time routing option by GA, focusing on a flight from MUC to JFK.
www.geosci-model-dev.net/9/3363/2016/Geosci.Model Dev., 9, 3363-3392, 2016 ues were similar and indicated a similar trend: an increased flight time for westbound flights on the trans-Atlantic region in winter.The mean fuel consumption values simulated by AirTraf were comparable to the reference values of BADA corresponding to low and nominal weights.The mean EINO x values were in the same range as the reference values of earlier studies.Finally, obtained maximum and minimum aircraft weights were compared to the three structural weight limits and one specified weight limit of the A330-301 aircraft.Almost all the values satisfied the four weight limits and only 15 of 515 flights exceeded the maximum take-off weight.Thus, AirTraf comprises a sufficiently good fuel use model.
The fundamental framework of AirTraf has been developed to perform fairly realistic air traffic simulations.Air-Traf 1.0 is ready for more complex routing tasks.Objective functions corresponding to other routing options will be integrated soon, and AirTraf will be coupled with various submodels of EMAC to evaluate air traffic climate impacts.

Code availability
AirTraf is published for the first time as a submodel of the Modular Earth Submodel System (MESSy).MESSy is continuously developed and used by a consortium of institutions.The usage of MESSy and access to the source code is licensed to all affiliates of institutions that are members of the MESSy Consortium.Institutions can become a member of the MESSy Consortium by signing the MESSy Memorandum of Understanding.More information can be found on the MESSy Consortium website (http:// www.messy-interface.org).The version presented here corresponds to AirTraf 1.0.Some improvements will be performed and AirTraf 1.0 will be updated for the latest version of the code.The status information for AirTraf including the license conditions is available on the website.
The Supplement related to this article is available online at doi:10.5194/gmd-9-3363-2016-supplement.

Figure 1
Figure1.Flowchart of EMAC/AirTraf.MESSy as part of EMAC provides interfaces (yellow) to couple various submodels for data exchange, run control and data input/output.Air traffic data and AirTraf parameters are input in the initialization phase (messy_initialize, dark blue).AirTraf includes the flying process in messy_global_end (dashed box, light blue), which comprises four main computation procedures (bold-black boxes).The detailed procedures are described in Sect.2.4 and are illustrated in Fig.3.AirTraf is linked to three modules: the aircraft routing module (light green), the flight trajectory optimization module (dark green), and the fuel/emissions calculation module (light orange).Resulting flight trajectories and global fields are calculated for output (rose red).Various submodels of EMAC can be linked to evaluate climate impacts on the basis of the output.

Figure 2 .
Figure2.Decomposition of global flight plans in a parallel environment of EMAC/AirTraf.A 1-day flight plan is distributed among many processing elements (PEs) in messy_init_memory (blue).A whole trajectory of an airport pair is handled by the same PE in the time loop of EMAC (messy_global_end, light blue).Finally, results are gathered from all the PEs for output (rose red).

Figure 3 .
Figure 3. Illustration of the flying process of AirTraf (dashed box in Fig. 1, light blue).(a) Flight trajectory calculation.(b) Fuel/emissions calculation.(c) Aircraft position calculation.(d) Gathering global emissions; the fraction of NO x,i corresponding to the flight segment i is mapped onto the nearest grid point (closed circle) relative to the (i + 1)th waypoint (open circle).ETO: estimated time over; F cr : fuel flow of an aircraft; m: aircraft weight; t: time step index of EMAC.The detailed calculation procedures are described in Sect.2.4.

Figure 4 .
Figure 4. Five representative routes for the great circle benchmark test.The details of locations are listed in Table3.

Figure 5 .
Figure 5.Comparison of the flight distance for the five representative routes.•: great circle distance calculated by Eq. (22); •: great circle distance calculated by Eq. (23).

Figure 6 .
Figure 6.Geometry definition of flight trajectory in the vertical cross section (top) and projection on the Earth (bottom).The bold solid line indicates a trajectory from MUC to JFK. •: control points determined by design variables x = (x 1 , x 2 , • • •, x 11 ) T .The lower/upper bounds of the 11 design variables are shown in Table6.Bottom: the dashed boxes show rectangular domains of three control points.: central points of the domains are calculated on the great circle (thin solid line), which divide the λ airport into four equal parts.Top: the dashed lines show the lower/upper variable bounds in altitude."FL290" stands for a flight level at 29 000 ft. Longitude coordinates for x 7 , x 8 , • • •, x 11 are pre-calculated; the coordinates divide the λ airport into six equal parts.

Figure 7 .
Figure 7. Optimal solutions varying with the population size n p and the number of generations n g .f means the difference in flight time between the optimal solution f and the true-optimal solution f true (= 25 994.0 s).The f (in %) is calculated as ( f/f true ) × 100.The flight time of the best solution is f best = 25 996.6 s (for n p = 100 and n g = 100, f < 3.0 s (less than 0.01 %)).

Figure 8 .
Figure8.Ten-thousand explored trajectories (solid line, black) from MUC to JFK in the vertical cross section (top) and projection on the Earth (bottom).The population size n p is 100 and the number of generations n g is 100.The best solution (red line) overlaps with the true-optimal solution (dashed line, black), i.e., the great circle at FL290.The flight time of the best solution is 25 996.6 s, while that of the true-optimal solution is 25 994.0 s.

Figure 9 .
Figure 9. Trajectories for the best solution (red line) and the trueoptimal solution (dashed line, black).This shows the enlarged drawing of Fig. 8 (top).The maximum difference in altitude is 0.83 m.

Figure 10 .
Figure 10.Flight time vs. number of function evaluations (= n p × n g ), including the enlarged drawing in the early 1000 evaluations.The population size n p is 100 and the number of generations n g is 100.f means the difference in flight time between the solution f and the true-optimal solution f true (= 25 994.0 s).The f (in %) is calculated as ( f/f true ) × 100.The solution shown as red line corresponds to the best solution in Figs.7 to 9. TableS1summarizes the 10 optimal solutions in detail.

Figure 11 .
Figure 11.Variation of the mean value of the difference in flight time between the true-optimal solution (f true = 25 994.0 s) and the optimal solution f (a), and the standard deviations of f (s f , b) are shown varying with the population size n p and the number of generations n g .The variation was calculated by 10 independent GA simulations from different initial populations for each combination of n p and n g : in total, 1000 independent simulations.On the f and time, UTC Mean flight altitude h, m (in ft)Flight time FT, s Time-optimal GC FL290 GC FL330 GC FL370

Figure 13 .
Figure 13.Ten-thousand explored trajectories (black lines) between MSP and AMS in the vertical cross section (top) and projection on the Earth (bottom), including the time-optimal flight trajectories (red and blue lines).(a) The eastbound flight from MSP to AMS.(b) The westbound flight from AMS to MSP.

Figure 14 .
Figure 14.Trajectories for the time-optimal (red and blue lines) and great circle cases (black lines) between MSP and AMS.The contours show the zonal wind speed (u in ms −1 ); arrows (black) show the wind speed ( u 2 + v 2 ) and direction.(a) The eastbound flight from MSP to AMS with the wind field at h = 8839 m at 21:35:00 UTC.(b) The westbound flight from AMS to MSP with the wind field at h = 10 002 m at 12:50:00 UTC.

Figure 15 .
Figure 15.Altitude distributions of the true air speed V TAS in ms −1 (a, b) and the tail wind indicator V ground /V TAS (c, d) along the time-optimal flight trajectories (black line) between MSP and AMS.Note that (V ground /V TAS ) ≥ 1.0 means tail winds (TW, red), while (V ground /V TAS ) < 1.0 means head winds (HW, blue) in the flight direction.The contours were obtained at the departure time: 21:35:00 UTC (eastbound, a, c); 12:50:00 UTC (westbound, b, d).

Figure 17 .
Figure 17.Obtained flight trajectories from 1-day AirTraf simulations corresponding to the time-optimal case including altitude changes in [FL290, FL410] (a, b) and the great circle cases at FL290, FL330, FL370 and FL410 (c, d).For each figure, the trajectories in the vertical cross section (top) and projection on the Earth (bottom).The 1-day flights comprise 52 eastbound (red lines) and 51 westbound flights (blue lines).

Figure 18 .
Figure18.Values of the true air speed V TAS (a) and the tail wind indicator V ground /V TAS (b) at waypoints for the time-optimal and great circle flights.Linear fits of the time-optimal (solid line, red (eastbound) and blue (westbound)) and great circle cases (dashed line, red (eastbound) and blue (westbound)) are included.V TAS of the international standard atmosphere (ISA) is given in (a) (solid line, black) provided by the BADA atmosphere table(Eurocontrol,  2010).

Figure 19 .
Figure 19.Mean fuel consumption (in kg(fuel) min −1 ) vs. altitude for the time-optimal and great circle flights.: mean value of all 103 flights; these values are shown in column 4 of Table10.

Figure 20 .
Figure 20.Comparison of the flight time for individual flights.A symbol indicates the value for one airport pair, corresponding to the time-optimal and great circle flights.If the value for the timeoptimal flight is the same as that of the great circle flight, the symbol lies on the 1 : 1 solid line.

Figure 21 .
Figure 21.Comparison of aircraft weights with structural weight limits (MTOW and MLW) and one specified weight limit (MLOW).The aircraft weights of the 103 flights for the time-optimal and great circle cases are plotted.•: aircraft weight at the last waypoint (m n wp ).•: aircraft weight at the first waypoint (m 1 ).The description of the limits is shown in Table14.

Table 1 .
Primary data of Airbus A330-301 aircraft and constant parameters used in AirTraf simulations.

Table 3 .
Information for the five representative routes of the great circle benchmark test.

Table 4 .
Great circle distance (d) of the five representative routes calculated with different calculation methods.Column 2 (d eq22 ) corresponds to the result calculated by Eq. (22); column 3 (d eq23 ) corresponds to the result calculated by Eq. (23) with n wp = 100; and column 4 (d MTS ) and d eq22,MTS =

Table 5 .
Calculation conditions for the benchmark test on flight trajectory optimizations.

Table 6 .
Lower/upper bounds of the 11 design variables.

Table 11 .
Flight time, fuel use, and NO x and H 2 O emissions for the time-optimal and great circle cases obtained from 1-day AirTraf simulations.Eastbound: sum of 52 eastbound flights; Westbound: sum of 51 westbound flights; and Total: sum of 103 flights.Changes (in %) relative to the time-optimal case are given in parentheses.