# Oceananigans.jl: A model that achieves breakthrough resolution, memory and
energy efficiency in global ocean simulations

Simone Silvestri  Massachusetts Institute of Technology, Cambridge, MA, USA
Gregory Wagner  Massachusetts Institute of Technology, Cambridge, MA, USA
Christopher Hill  Massachusetts Institute of Technology, Cambridge, MA, USA
Matin Raayai Ardakani  Northeastern University, Boston, MA, USA  Johannes
Blaschke  Lawrence Berkeley National Laboratory, Berkeley, CA, USA  Jean-
Michel Campin  Massachusetts Institute of Technology, Cambridge, MA, USA
Valentin Churavy  Massachusetts Institute of Technology, Cambridge, MA, USA
Navid Constantinou  Australian National University, Canberra, ACT, Australia
Alan Edelman  Massachusetts Institute of Technology, Cambridge, MA, USA  John
Marshall  Massachusetts Institute of Technology, Cambridge, MA, USA  Ali
Ramadhan  Massachusetts Institute of Technology, Cambridge, MA, USA  Andre
Souza  Massachusetts Institute of Technology, Cambridge, MA, USA  Raffaele
Ferrari  Massachusetts Institute of Technology, Cambridge, MA, USA

(March 2023)

###### Abstract

Climate models must simulate hundreds of future scenarios for hundreds of
years at coarse resolutions, and a handful of high resolution decadal
simulations to resolve localized extreme events. Using Oceananigans.jl,
written from scratch in Julia, we report several achievements: First, a global
ocean simulation with breakthrough horizontal resolution — 488m — reaching 15
simulated days per day (0.04 simulated years per day; SYPD). Second,
Oceananigans simulates the global ocean at 488m with breakthrough memory
efficiency on just 768 Nvidia A100 GPUs, a fraction of the resources available
on current and upcoming exascale supercomputers. Third, and arguably most
significant for climate modeling, Oceananigans achieves breakthrough energy
efficiency reaching 0.95 SYPD at 1.7 km on 576 A100s and 9.9 SYPD at 10 km on
68 A100s — the latter representing the highest horizontal resolutions employed
by current IPCC-class ocean models. Routine climate simulations with 10 km
ocean components are within reach.

##  1 Justification

Oceananigans.jl — a new ocean model written from scratch in Julia — achieves
ocean simulations with breakthrough resolution, memory and energy efficiency,
realizing 0.041 simulated years per day (SYPD) at 488 m on 768 Nvidia A100s,
0.95 SYPD at 1 km on 576 A100s, and 9.9 SYPD at 10 km on 68 A100s.

##  2 Performance Attributes

Categories | Scalability, time-to-solution, energy-to-solution.
---|---
Type of method | Fully explicit with sub-cycling.
Results basis | Whole application excluding I/O.
Numerical precision | Both 64- and 32-bit cases measured.
System scale | Results measured on full-scale systems.
Measurement mechanism | Timers, memory used and energy used.

##  3 Overview of the Problem

Climate models are essential for predicting where, when, and how climate
change threatens Earth’s ecosystems and human civilization. But current
climate models, which capture only the broadest aspects of global warming,
fall far short of providing the needed accuracy and granularity required to
design and implement costly adaptation and mitigation strategies [14].
Significant reduction of the uncertainty of climate predictions is potentially
worth trillions of dollars [20].

Climate models simulate the three-dimensional fluid dynamics, thermodynamics,
chemistry, and biology of the atmosphere, ocean, and land to predict the
hydrological cycle, carbon cycle and the net energy imbalance of the Earth
system. While typical climate models use coarse resolutions of 25-100 km to
simulate the numerous climate scenarios required by the Intergovernmental
Panel on Climate Change (IPCC) [25], a handful of state-of-the-art climate
simulations have been performed at higher resolutions of O(10 km) at
astronomical expense. At either resolution there are many processes, such as
clouds and ocean turbulence, that cannot be explicitly simulated and are
instead approximated by empirical formulae called parameterizations. Biases
due to inadequate parameterizations dominate the uncertainty of climate
predictions over the next few decades [34, 23].

The prevailing strategy to reduce climate model uncertainty is to refine model
resolution as much as possible [34]. For example, at horizontal resolutions of
1 km a substantial fraction of atmospheric convection and ocean turbulence are
explicitly modeled by Newton’s laws of motion, greatly reducing the impact of
parameterizations [34]. High-resolution climate modeling is further required
to make predictions for specific regions, providing information for local
decision makers on adaptation and mitigation [14].

Yet the “resolution strategy” is fundamentally limited: even at 1 km
resolution many climate-relevant physical processes remain unresolved [49].
Worse, processes such as sea ice dynamics, biology, or cloud-aerosol
interaction will never be resolved because accurate macroscopic laws do not
exist. Absent theoretical breakthroughs, such “irreducible” uncertainties can
be addressed only by leveraging Earth system observations through advances in
data assimilation and machine learning [41]. Data-driven optimization of
climate models requires ensembles of climate predictions, rather than single
predictions at the highest affordable resolution. Ensembles of simulations are
also required to explore emission scenarios and to estimate the impact of
initial condition uncertainty and internal variability.

Consequently, reducing the uncertainty of climate predictions demands not just
higher resolution, but more efficient resource utilization to enable hundreds
to thousands of relatively high-resolution simulations. As an example, we
consider the computational requirements to enable 100-simulation ensembles
using all 37,888 AMD MI250 GPUs of the Frontier exascale supercomputer:
completing an ensemble of 300-year simulations (200 years of spin-up + 100
years of prediction) within one month of wall clock time requires a climate
model that can achieve 10 simulated years per day (SYPD) using 378 GPUs, or
1/100th of Frontier’s resources. Disruptive progress on climate modeling
requires not just scalable performance for a single, high-resolution
simulation, but advances in efficiency to meet this ensemble-based “10 per
100th” benchmark [40].

Our submission uses the ocean component of a new climate model being developed
by the Climate Modeling Alliance [44]. The ocean contributes key uncertainty
to climate predictions due to its prominent role in the Earth system’s heat
and carbon cycles. At 10 km resolutions, where ocean model uncertainties are
significantly reduced, the ocean often is the most expensive climate model
component [19]. This calls for a step change in ocean model performance.

##  4 Current State of the Art

We are aware of only three global ocean simulations that have achieved
resolutions finer than 5 kilometers — all at tremendous computational expense.
In 2014, MITgcm [29] was used to perform the one year, tidal-forced ice-ocean
simulation “LLC4320” [45], which exhibits 2.2 km horizontal resolution with 90
vertical levels. LLC4320 achieved 0.047 simulated years per day (SYPD) using
70,000 cores of the NASA Pleiades system.

Figure 1: Simulated years computed by a megawatt-hour of energy (SWPMWh)
versus number of grid points for state-of-the-art atmosphere and ocean models.
Stars show the performance of our ocean model in a realistic and “aqua planet”
(AP) setup.

FIO-COM32 [50] ran at $\sim$2.5 km (1/32nd degree) horizontal resolution with
90 vertical levels for 3.5 years. [48] ported LICOM3 to GPUs to realize 0.51
SYPD at 1/20th horizontal resolution with 60 vertical levels using 384 MI50
AMD GPUs, and further managed to scale to 26200 MI50s with strong scaling
efficiency of 8%.

The largest ocean simulations used in current IPCC-class climate models, which
typically require faster-time-to-solution to support longer simulations, have
horizontal resolutions of roughly 10 km. [8] describes output from four
60-year ocean simulations following the OMIP-2 protocol with 8 km (1/12th
degree), 10 km, and two with 11 km (1/10th degree). [11] report a 110-year
simulation at 10 km (1/10th degree) horizontal resolution, the longest high
resolution OMIP-2-style run. Some of the highest resolution climate models are
the iHESP CESM-based model with 25km-10km atmosphere-ocean resolution [51],
achieving 3.4 SYPD, and the 50km-10km HadGEM3-GC3.1 submission to HighResMIP
[18, 37], achieving 0.4 SYPD.

At 3.4 SYPD, the iHESP CESM achieves sufficient time-to-solution for hundreds
to thousands of years of simulated climate. But such a simulation is purchased
for a high price, requiring the 40% of the Sunway TaihuLight supercomputer
[51] and 4 million cores consuming 6 MW for hundreds of days of wall clock
time. Enabling the large ensembles of high-resolution simulations needed to
improve climate prediction requires both performance at scale as well as
efficient resource utilization.

Figure 1 plots simulated years per mega-watt-hour (SYPMWh) against resolution
for state-of-the-art ocean models. The SYPMWh metric encodes the efficiency
requirement needed to make progress on climate uncertainty with next-
generation climate models: in particular, we require both higher-resolution
models (moving rightwards in figure 1) as well as more efficient models
(moving upwards in figure 1). For completeness, we report SYPMWh also for two
GPU-based models: Veros [22], an ocean model, and COSMO [15], an atmospheric
model. The present nomination is shown with stars from whence we see
significant performance gains compared to the existing state-of-the art.

##  5 Innovations

Our achievement is three-fold: first, using new software written in the Julia
programming language called Oceananigans.jl [35], we report a near-global
ocean simulation with highest-ever horizontal resolution (488 m) reaching 15
simulated days per day (0.04 SYPD). Second, Oceananigans performs this
simulation with breakthrough memory efficiency on just 768 NVidia A100 GPUs,
and thus a fraction of the available resources on current and upcoming
exascale supercomputers. Third, and arguably most important, Oceananigans
achieves breakthrough energy efficiency, simulating the global ocean at 0.95
SYPD with 1.7 km resolution on 576 A100s, and at 10 km — the highest
horizontal resolution employed by an IPCC-class ocean model — achieving 9.9
SYPD on 68 Nvidia A100s. This final milestone proves the feasibility of
routine climate simulations with 10 km ocean components, a crucial resolution
threshold at which ocean macroturbulence (the most energetic ocean motions
with scales between 10–100 km) is fully resolved.

We attribute these achievements first and foremost to a high-risk, high-reward
strategy to develop a new ocean model from scratch in Julia with a specific
focus on GPU performance and memory efficiency. Additional crucial ingredients
include advances in numerical methods for finite volume fluid dynamics on the
sphere and a novel optimization for simulating ocean free surface dynamics
that achieves unprecedented GPU scalability.

###  5.1 Starting from scratch with Julia

Oceananigans.jl is an open-source library for ocean-flavored fluid dynamics
written from scratch in Julia [7]. Julia is a dynamic high-level programming
language that leverages Just-In-Time (JIT) compilation and LLVM [24] to
achieve performance competitive with traditional HPC languages like C or
Fortran. Julia has gathered interest as potential language for HPC [17, 10,
16, 21, 27] and provides easy integration with MPI [47, 38]. Most of
Oceananigans.jl software is hardware-agnostic through the Julia package
KernelAbstractions.jl [10], which enables performance portability targeting
CPUs and different GPU vendors using the JuliaGPU [6, 5] software stack,
similar to the capabilities provided by Kokkos [9], OCCA [30], and HIP [1].

To our knowledge, Oceananigans is the first ocean model written from scratch
for GPUs, rather than ported from existing CPU code. Starting from scratch and
using the Julia programming language allowed us to rethink the typical
patterns used in ocean and atmosphere dynamical cores. In particular, we
developed a system of composable atomic operators that leverages Julia’s
functional programming paradigm and effective inlining capabilities to
recursively construct large expression trees for calculus on staggered finite
volume grids. Using this composable operator system, we fuse the entire
tendency computation for each prognostic variable into a single compute-heavy
kernel, each of which depends on only two intermediate diagnostic variables
representing hydrostatic pressure and vertical diffusivity (which is treated
implicitly using a predictor-corrector method).

Such a high degree of abstraction yields a number of innovations: first,
kernel fusion maximizes efficiency on GPUs. Second, almost all intermediate
quantities are computed on-the-fly, so that Oceananigans is extremely memory
efficient and can perform global ocean simulations at resolutions up to 1/4th
degree on a single Nvidia V100. Finally, because all compute-heavy kernels
rely on a single “tendency kernel function” applied at each grid index
$i,j,k$, we can easily optimize performance by rapidly prototyping techniques
to overlap computation and communication. The sparsity of kernels per time-
step and small number of temporary variables mean that Oceananigans’
algorithmic structure is markedly different from current ocean models, which
typically allocate 10 to 100 times the minimum necessary memory [4] and
distribute computations across many small kernels [51]. We argue these
algorithmic differences are a major factor in Oceananigans’ energy-efficiency
and time-to-solution on GPU systems.

###  5.2 New numerical methods for finite volume fluid dynamics on the sphere

Our results use Oceananigans.HydrostaticFreeSurfaceModel, which solves the
hydrostatic Boussinesq equations in a finite volume framework on staggered
C-grids [3]. Oceananigans’ hydrostatic model employs an implicit-explicit
second-order Adams-Bashforth time stepping scheme. Vertically-implicit
diffusion is implemented with a backward Euler time-discretization and
tridiagonal solver.

A major innovation is a new adaptive-order scheme based on weighted
essentially non-oscillatory (WENO) reconstructions [42] for advecting momentum
and tracers on curvilinear finite-volume grids [43]. This new scheme
automatically adapts to changing spatial resolution and permits stable, high-
fidelity simulations of ocean turbulence without explicit dissipation or
hyper-dissipation. This innovation reduces setup time when changing or
increasing resolution while guaranteeing high-fidelity solutions that exhibit
the minimum necessary dissipation of sharp, near-grid scale features.

Figure 2: Left: time-stepping sequence. Right: different domains over which 2D
fast and 3D slow mode updates take place (here assuming 1 barotropic substep
per baroclinic step – halo region of size 1 – and second-order methods – outer
region of size 1)

###  5.3 Optimization of ocean free surface dynamics for unprecedented GPU
scalability

In hydrostatic ocean models with a free surface, the vertically-averaged, two-
dimensional “barotropic mode” has dynamics orders of magnitude faster than the
three-dimensional “baroclinic” component, and must be treated by a special
“barotropic solver”. Due to communication overhead, barotropic solvers in
current ocean models — whether implicit or explicit — are a major bottleneck
that accounts for between 40% [22] to 60% [48, 36] of the cost of a typical
IPCC-class ocean simulations.

Oceananigans’ excellent scalability is enabled by an innovative optimization
of the parallel barotropic solver. An increase in computation is traded in for
decreased communication latency by leveraging the two-dimensionality of the
barotropic problem. Our new barotropic solver is based on explicit subcycling
of the barotropic mode. Increasing the width of the barotropic halo to equal
the number of explicit subcycles (typically between 10–30) greatly decreases
the frequency of communication. As a result, communication is required once
per time-step rather than every subcycle, reducing the frequency of
communication by a factor of 10 to 30. The cost of the barotropic solver is
therefore less than 10% of the total cost of a time step. Due to the sparsity
of communication enabled by our novel barotropic solver, all communication
operations can be overlapped with computational workloads as sketched in
figure 2.

##  6 How performance was measured

The Oceananigans model performance is estimated for two near-global ocean
simulations with different domains: a realistic (R) domain and an aqua planet
(AP) domain. Both domains span the entire longitudinal extent of the sphere
and cover a latitude range of $75^{\circ}$S to $75^{\circ}$N.

The Realistic domain has realistic bathymetry and is forced by realistic
surface momentum, heat, and salinity fluxes derived from the ECCO2 state
estimate[31] at three resolutions:

  * •

OceananigansR12 with 1/12th degree horizontal resolution ($\sim$7 km) and 48
vertical levels

  * •

OceananigansR24 with 1/24th degree horizontal resolution ($\sim$3.4 km) and
100 vertical levels

  * •

OceananigansR48 with 1/48th degree horizontal resolution ($\sim$1.7 km) and
100 vertical levels

Figure 3 shows surface vertical vorticity after one year integration of
OceananigansR12 and OceananigansR48 over the global ocean and also for
selected regions to show further detail. Both OceananigansR48 and
OceananigansR12 exhibit macroscale turbulent ocean features that are currently
unresolved by most IPCC-class models. The OceananigansR48 solution exhibits
fronts, filaments, and other “submesoscale” vorticity features realized only a
handful of times in global simulations.

The idealized OceananigansAP suite of simulations [13], which has idealized
bathymetry and surface forcing that does not require interpolation to
different resolutions, is used for weak scaling experiments. All
OceananigansAP experiments have 100 vertical levels and two latitudinal ridges
that divide the world ocean into two basins. We vary the horizontal resolution
of OceananigansAP from 1/6th of a degree ($\sim$14 km) to 1/196th of a degree
($\sim$488 m).

None of our simulations require explicit horizontal diffusion of momentum or
tracers owing to the adaptive WENO advection scheme described in section 5.2.
All simulations use a Richardson-number-based parameterization for vertical
mixing due to unresolved shear and convective turbulence at 1–100 m scales.

To assess the time-to-solution for each experiment in simulated years per day
(SYPD), we measure average wall clock time per time-step. Wall clock time is
sampled through NVIDIA’s Nsight System and recorded by NVIDIA Tool Extension
Library via the NVTX.jl Julia package.

To assess the efficiency of each solution in simulated years per mega-watt-
hour (SYPMWh), we combine SYPD with an estimate of the mean power draw over
the duration of an experiment. On MIT Satori [2], which has 256 Nvidia V100s,
we have access to precise, billing-grade power metering. For all simulations
with Nvidia A100s we estimate power consumption $P$ with

$P=250D+300N\,\text{Watts}\,,$ (1)

where $D$ is the number of A100s and $N$ is the number of nodes.

We further note that power estimates are provided by LICOM3 and COSMO, but not
for LLC4320 or Veros. To estimate the power consumption of LLC4320, we assume
that each of the 1000 dual CPU nodes draws 500W. We estimate the power
consumption of iHESP CESM [51] and HadGCM3 [37] as a percentage of the peak
power consumption of their respective clusters. We use equation (1) to
estimate Veros’ power consumption on 1 node with 16 A100s.

Figure 3:  Vertical vorticity as simulated by Oceananigans12 (top left) and
Oceananigans48 (bottom left) after a one year integration on September 1st. To
the right, insets zoom on particularly energetic current systems: the Aghulas
and the East Australian Currents. While major ocean currents with widths of
10-100 km are resolved in both simulations, the sharp density fronts and
associated currents that develop at the ocean surface in winter at scales
between 1-10 km (the ocean weather) are only resolved by Oceananigans48. On
September 1 — spring in the southern hemisphere, fall in the northern
hemisphere — such sharp frontal features populate the southern ocean but are
suppressed in the north.

##  7 Performance Results

We report both scaling results via time-to-solution in SYPD and efficiency
results via energy-to-solution in SYPMWh.

###  7.1 Scaling Results

Figure 4: Strong scaling tests for the realistic setups OceananigansR12
($1/12^{\circ}$), OceananigansR24 ($1/24^{\circ}$), and OceananigansR48
($1/48^{\circ}$). The left plot reports simulated years per wall clock day
(SYPD) while the right plot wall clock milliseconds per time steps. All
results are averaged over 1500 time steps.

Realistic ocean simulations (Satori and Engaging clusters). We report strong
scaling tests using the realistic global setup shown in figure 3 on two
clusters: (i) the MIT Satori cluster [2], a high-performance Power 9 system
composed of 64 Power 9 nodes hosting four Nvidia V100 GPUs with 32GBs memory
each, and (ii) the Engaging MIT cluster, using 8 nodes that host 4 NVlinked
A100s with 80GBs memory each. The resulting wall clock time per time step,
averaged over 1500 time steps, is presented in Figure 4 for both single
precision (FP32) and double precision (FP64) computations. On a single node,
OceananigansR12 attains 0.9 SYPD in double precision and 1.4 SYPD in single
precision, with a wall clock time per time step ranging from 330 to 550 ms.
When increasing the number of nodes up to 16 (64 GPUs), the communication
overhead increases, resulting in 12.4 SYPD in single precision and 7.75 SYPD
in double precision. We measure a strong scaling efficiency of 52% in single
precision and 55% in double precision over 64 GPUs, because the computational
workload (40 ms wall clock time per time-step) eventually becomes too short to
completely mask the communication overhead.

For higher-resolution ocean weather-permitting simulations, the scaling is
almost ideal across the range we investigate. For OceananigansR24 (FP64-V100)
and OceananigansR48 (FP32-V100), we measure larger than ideal scaling. This
counter-intuitive result is a product of a load balance improvement as the
number of GPUs increases. In summary, we attain 1.94 SYPD on 120 V100 GPUs
with a kilometer-scale resolution (OceananigansR24) and 0.33 SYPD with an
ocean weather-resolving simulation (OceananigansR48). Finally, we have tested
the OceananigansR48 setup on 144 Perlmutter nodes (576 A100 GPUs), reaching
the 0.95 SYPD. This is the first instance of a kilometer-scale ocean achieving
$\sim$1 SYPD. We have also tested the OceananigansR12 setup on 17 nodes
obtaining 9.9 SYPD (see fig. 5).

Figure 5: Weak scaling tests performed in double precision with the
OceananigansAP setup. Each GPU has a grid equivalent to a global $1/6^{\circ}$
and 100 vertical layers. The weak scaling is performed up to a horizontal
resolution of 1/168th of a degree ($\sim$488 m resolution) where we achieve 15
simulated days per wall clock day (1 year in roughly 25 days). The star marks
the performance of OceananigansR48 (figure 3) on 144 Perlmutter GPU nodes. All
results are averaged over 500 time steps.

Aqua-planet simulation (Perlmutter cluster). We report weak scaling tests on
the NERSC supercomputer (Perlmutter). Perlmutter is a HPE (Hewlett Packard
Enterprise) Cray EX supercomputer that hosts four A100 GPUs with 40GB per
node, linked through a NVLink3 interconnect. All weak scaling tests are
performed using the OceananigansAP setup on double precision. We allocate two
different horizontal resolutions (1/12 and 1/6 of a degree), progressively
increasing them with the number of GPUs while maintaining 100 vertical levels.
As shown in figure 5, we obtain 100% weak scaling efficiency for the whole
investigated range (1 to 196 nodes – 4 to 768 A100s).

###  7.2 Energy efficiency

In table 1 we summarize the energy metrics for our computations as well as the
other investigated models. Figure 1 is derived from the data outlined in this
table. HadGEM3 and iHRES entries are estimated by including the whole coupled
climate model (atmosphere and ocean). Unavailable data is marked with $-$. Our
Oceananigans simulations are the highest in each of their columns. This
reflects our attention to memory and energy efficiency.

Model | Time step | Grid size | CPU/GPU | SYPD | wtime/tstep | Power est. | In fig 1
---|---|---|---|---|---|---|---
HadGEM3${}_{\text{FP64}}$ (Climate) [37] | - | $\sim 7.22\times 10^{8}$ | 9396 (Cray XC40) | 0.4 | - | 141KW | $\usym{1F5F8}$
iHRES${}_{\text{FP64}}$ (Climate) [51] | - | $\sim 6\times 10^{8}$ | Sunway TaihuLight | 3.7 | - | 6500KW | $\usym{1F5F8}$
LLC4320${}_{\text{FP64}}$ (Ocean) | 25 s | $8.7\times 10^{10}$ | 2000 (Intel) | 0.041 | 1.6 | 500KW | $\usym{1F5F8}$
Veros${}_{\text{FP64}}$ (Ocean) [22] | 180 s | $3.5\times 10^{8}$ | 16 (A100) | 0.8 | 0.62 | 5.2KW | $\usym{1F5F8}$
Veros${}_{\text{FP32}}$ (Ocean) [22] | 180 s | $3.5\times 10^{8}$ | 16 (A100) | 1.3 | 0.38 | 5.2KW |
LICOM3${}_{\text{FP64}}$ (Ocean) [48] | 60 s | $1.5\times 10^{9}$ | 384 (MI50) | 0.51 | 0.32 | 92KW | $\usym{1F5F8}$
LICOM3${}_{\text{FP64}}$ (Ocean) [48] | 60 s | $1.5\times 10^{9}$ | 26200 (MI50) | 2.72 | 0.06 | 6300KW |
COSMO${}_{\text{FP64}}$ (Atmos) [15] | 6 s | $3.46\times 10^{10}$ | 4888 (P100) | 0.043 | 0.4 | 1000KW | $\usym{1F5F8}$
OceananigansR12${}_{\text{FP32}}$ (Ocean) | 180 s | $3.7\times 10^{8}$ | 4 (V100) | 1.5 | 0.33 | 1.2KW | $\usym{1F5F8}$
OceananigansR12${}_{\text{FP32}}$ (Ocean) | 180 s | $3.7\times 10^{8}$ | 64 (V100) | 12.4 | 0.04 | 18KW |
OceananigansR12${}_{\text{FP64}}$ (Ocean) | 180 s | $3.7\times 10^{8}$ | 68 (A100) | 9.9 | 0.05 | 22 | $\usym{1F5F8}$
OceananigansR48${}_{\text{FP32}}$ (Ocean) | 45 s | $1.24\times 10^{10}$ | 120 (V100) | 0.33 | 0.37 | 36KW |
OceananigansR48${}_{\text{FP64}}$ (Ocean) | 45 s | $1.24\times 10^{10}$ | 576 (A100) | 1.0 | 0.13 | 187KW | $\usym{1F5F8}$
OceananigansR48${}_{\text{FP64}}$ (Ocean) | 45 s | $1.24\times 10^{10}$ | 32 (A100) | 0.063 | 1.9 | 10.4KW |
OceananigansAP${}_{\text{FP64}}$ (Ocean) | 11 s | $2.1\times 10^{11}$ | 768 (A100) | 0.063 | 0.81 | 252KW | $\usym{1F5F8}$

Table 1:  Performance details of state-of-the are climate, ocean, and
atmosphere models. Larger grid sizes correspond to finer spatial resolution.
Computations belonging to this submission are shown in bold.

##  8 Implications

By developing a new model from scratch specifically for GPUs, and wielding a
handful of key ocean-model-specific innovations, Oceananigans achieves 9.9
SYPD at 10 km resolution using less than 1% of the resources of current state
of the art supercomputers. This achievement means that most climate model runs
submitted to IPCC will be able use 10 km ocean models — precipitating a step
change in the accuracy of climate prediction.

At scales between 10–100 km, macroscale ocean turbulence exerts a key control
on ocean carbon and heat uptake. However, attempts to accurately parameterize
this key process in coarse resolution models have frustrated generations of
oceanographers. The inadequacies of macroscale parameterizations are
associated with major biases and uncertainty in climate predictions [33, 39].
At resolutions of 10 km, the need for macroscale turbulence parameterization
is eliminated, and ocean simulations capture key ocean features such as sharp
sea surface temperature gradients supporting the formation of marine stratus
clouds above narrow eastern boundary currents like the California and Benguela
Current [28], and changes in the meridional overturning circulation due to the
effect of Antarctic meltwater on deep convection in austral winter [26].

Additionally, by achieving 0.95 SYPD at 1.7 km resolution, we pave the way for
decadal ocean simulations of the ocean “submesoscale” — the ocean analogue to
atmospheric weather — which exhibits hourly fluctuations, high spatial and
seasonal variability, and which exerts a strong control on ocean air-sea
fluxes, biological productivity and fish stocks [46]. The granularity and
accuracy provided by 1.7 km resolution is further required to plan local
mitigation strategies and predict local extreme events.

Third, the unparalleled speed of execution and memory efficiency of
Oceananigans allows global computations at never-before-seen sub-kilometer
resolutions. The capacity for ultra-high-resolution simulations aligns with
current advancements in resolution of ocean sampling platforms from satellites
[32, 12] to fleets of floats and drones. While this wealth of data is likely
to provide new insights and scientific knowledge about the nature of small
scale processes, global high-resolution ocean simulations will be needed to
explore their impact on global climate scales.

Finally, our results pave the way for marked increase in energy efficiency of
climate simulations. The very reason to develop climate models, as stated by
the Coupled Model Intercomparison Project (CMIP), for example, is to provide
the necessary information to effectively reduce emissions and mitigate the
effects of global warming — while, counterproductively, the carbon footprint
of climate simulations that contribute to CMIP increases rapidly.
Oceananigans’ achievements represent a milestone towards decreased energy
consumption by climate modeling efforts.

##  9 Acknowledgments

This research used resources of the National Energy Research Scientific
Computing Center (NERSC), a U.S. Department of Energy Office of Science User
Facility located at Lawrence Berkeley National Laboratory, operated under
Contract No. DE-AC02-05CH11231 using NERSC award DDR-ERCAP0025591. This work
is partly supported by the generosity of Eric and Wendy Schmidt by
recommendation of the Schmidt Futures program and by NSF grant AGS-1835576.
N.C.C. is supported by the Australian Research Council DECRA Fellowship
DE210100749.

## References

  * [1] HIP: C++ Heterogeneous-Compute Interface for Portability.  https://github.com/ROCm-Developer-Tools/HIP.
  * [2] MIT Satori User Documentation.  https://mit-satori.github.io/index.html.
  * [3] A. Arakawa and V.R. Lamb.  Computational design of the basic dynamical processes of the UCLA general circulation model.  General Circulation Models of the Atmosphere, 17, 1977.
  * [4] V. Balaji, E. Maisonnave, N. Zadeh, B.N. Lawrence, J. Biercamp, U. Fladrich, G. Aloisio, R. Benson, A. Caubel, J. Durachta, et al.  CPMIP: measurements of real computational performance of Earth system models in CMIP6.  Geoscientific Model Development, 10(1):19–34, 2017.
  * [5] T. Besard, V. Churavy, A. Edelman, and B. De Sutter.  Rapid software prototyping for heterogeneous and distributed platforms.  Advances in engineering software, 132:29–46, June 2019.
  * [6] T. Besard, C. Foket, and B. De Sutter.  Effective Extensible Programming: Unleashing Julia on GPUs.  IEEE Trans. Parallel Distrib. Syst., (4):827–841, December 2017\.
  * [7] J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah.  Julia: A fresh approach to numerical computing.  SIAM Review, 59(1):65–98, January 2017.
  * [8] E.P. Chassignet, S.G. Yeager, B. Fox-Kemper, A. Bozec, F. Castruccio, G. Danabasoglu, C. Horvat, W.M. Kim, N. Koldunov, Y. Li, et al.  Impact of horizontal resolution on global ocean–sea ice model simulations based on the experimental protocols of the Ocean Model Intercomparison Project phase 2 (OMIP-2).  Geoscientific Model Development, 13(9):4595–4637, 2020.
  * [9] T.R. Christian, D. Lebrun-Grandie, D. Arndt, J. Ciesko, V. Dang, N. Ellingwood, R. Gayatri, E. Harvey, D.S. Hollman, D. Ibanez, et al.  Kokkos 3: Programming model extensions for the exascale era.  IEEE Transactions on Parallel and Distributed Systems, 33(4):805–817, 2021.
  * [10] V. Churavy, D. Aluthge, J. Samaroo, A. Smirnov, J. Schloss, L. C Wilcox, S. Byrne, M. Waruszewski, A. Ramadhan, S. Schaub, N. C Constantinou, J. Bolewski, M. Ng, T. Besard, B. Arthur, C. Kawczynski, C. Hill, C. Rackauckas, J. Cook, J. Liu, M. Schanen, O. Schulz, P. Haraldsson, T. Arakaki, and T. Chor.  JuliaGPU/KernelAbstractions.jl: v0.9.1, March 2023.
  * [11] M. Ding, H. Liu, P. Lin, Y. Meng, W. Zheng, B. An, Y. Luan, Y. Yu, Z. Yu, Y. Li, J. Ma, J. Chen, and K. Chen.  A century-long eddy-resolving simulation of global oceanic large- and mesoscale state.  Sci. Data, 9:661, 2022.
  * [12] C. Donlon, B. Berruti, S. Mecklenberg, J. Nieke, H. Rebhan, U. Klein, A. Buongiorno, C. Mavrocordatos, J. Frerick, and B. Seitz.  The sentinel-3 mission: Overview and status.  In 2012 IEEE International Geoscience and Remote Sensing Symposium, pages 1711–1714. IEEE, 2012.
  * [13] D. Ferreira, J. Marshall, and J.-M. Campin.  Localization of deep water formation: Role of atmospheric moisture transport and geometrical constraints on ocean circulation.  J. Climate, 23:1456–1476, 2010.
  * [14] T. Fiedler, A.J. Pitman, K. Mackenzie, N. Wood, C. Jakob, and S.E. Perkins-Kirkpatrick.  Business risk and the emergence of climate analytics.  Nature Climate Change, 11:87–94, 2021.
  * [15] O. Fuhrer, T. Chadha, T. Hoefler, G. Kwasniewski, X. Lapillonne, D. Leutwyler, D. Lüth, C. Osuna, C. Schär, T.C. Schulthess, and H. Vogt.  Near-global climate simulation at 1 km resolution: establishing a performance baseline on 4888 GPUs with COSMO 5.0.  Geosci. Model Dev., 14:2781–2799, 2021.
  * [16] M. Giordano, M. Klöwer, and V. Churavy.  Productivity meets performance: Julia on A64FX.  In 2022 IEEE International Conference on Cluster Computing (CLUSTER), pages 549–555, September 2022.
  * [17] W. F. Godoy, P. Valero-Lara, E. T. Dettling, C. Trefftz, I. Jorquera, T. Sheehy, R. G. Miller, M. Gonzalez-Tallada, J. S Vetter, and V. Churavy.  Evaluating performance and portability of high-level programming models: Julia, Python/Numba, and kokkos on exascale nodes.  March 2023.
  * [18] J M. Gutiérrez and A.-M. Tréguier.  IPCC, 2021: Annex II: Models.  Climate Change 2021: The Physical Science Basis, page 2087–2138, 2021.
  * [19] R. Haarsma, M. Acosta R., Bakhshi, P.-A. Bretonnière, L.-P. Caron, M. Castrillo, S. Corti, P. Davini, E. Exarchou, F. Fabiano, U. Fladrich, R. Fuentes Franco, J. García-Serrano, J. von Hardenberg, T. Koenig, X. Levine, V.L. Meccia, T. van Noije, G. van den Oord, F.M. Palmeiro, M. Rodrigo, Y. Ruprich-Robert, P. Le Sager, E. Tourigny, S. Wang, M. van Weele, and K. Wyser.  HighResMIP versions of EC-Earth: EC-Earth3P and EC-Earth3P-HR – description, model computational performance and basic validation.  Geoscientific Model Development, 13(8):3507–3527, 2020.
  * [20] C. Hope.  The $10 trillion value of better information about the transient climate response.  Phil. Trans. R. Soc. A, 3:20140429, 2015.
  * [21] S. Hunold and S. Steiner.  Benchmarking Julia’s communication performance: Is Julia HPC ready or full HPC?  In 2020 IEEE/ACM Performance Modeling, Benchmarking and Simulation of High Performance Computer Systems (PMBS), pages 20–25. IEEE, November 2020.
  * [22] D. Häfner, R. Nuterman, and M. Jochum.  Fast, Cheap, and Turbulent—Global Ocean Modeling With GPU Acceleration in Python.  JAMES, 13(12):e2021MS002717, 2021.
  * [23] E.J. Kendo, N.M. Roberts, H.J. Fowler, M.J. Roberts, S.C. Chan, and C.A. Senior.  Heavier summer downpours with climate change revealed by weather forecast resolution model.  Nature Climate Change, 4(7):570–576, 2014.
  * [24] C Lattner and V Adve.  LLVM: A compilation framework for lifelong program analysis & transformation.  In International Symposium on Code Generation and Optimization, 2004\. CGO 2004. IEEE, 2004.
  * [25] H. Lee and J. Romero (eds.).  IPCC, 2023: Climate Change 2023: Synthesis Report. A Report of the Intergovernmental Panel on Climate Change. Contribution of Working Groups I, II and III to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change, 2023.
  * [26] Q. Li, M.H. England, A. McC Hogg, S.R. Rintoul, and A.K. Morrison.  Abyssal ocean overturning slowdown and warming driven by Antarctic meltwater.  Nature, 615(7954):841–847, 2023.
  * [27] W. C. Lin and S. McIntosh-Smith.  Comparing julia to performance portable parallel programming models for HPC.  In 2021 International Workshop on Performance Modeling, Benchmarking and Simulation of High Performance Computer Systems (PMBS), pages 94–105. ieeexplore.ieee.org, November 2021.
  * [28] J. Ma, S. Xu, and B. Wang.  Warm bias of sea surface temperature in eastern boundary current regions—a study of effects of horizontal resolution in CESM.  Ocean Dynamics, 69:939–954, 2019.
  * [29] J. Marshall, A. Adcroft, C. Hill, L. Perelman, and C. Heisey.  A finite-volume, incompressible Navier Stokes model for studies of the ocean on parallel computers.  Journal of Geophysical Research: Oceans, 102(C3):5753–5766, 1997\.
  * [30] D.S. Medina, A. St-Cyr, and T. Warburton.  OCCA: A unified approach to multi-threading languages.  arXiv preprint arXiv:1403.0968, 2014.
  * [31] Dimitris Menemenlis, Jean-Michel Campin, Patrick Heimbach, Chris Hill, Tong Lee, An Nguyen, Michael Schodlok, and Hong Zhang.  Ecco2: High resolution global ocean and sea ice data synthesis.  Mercator Ocean Quarterly Newsletter, 31(October):13–21, 2008.
  * [32] R. Morrow, L-L. Fu, F. Ardhuin, M. Benkiran, B. Chapron, E. Cosme, F. d’Ovidio, J.T. Farrar, S.T. Gille, G. Lapeyre, P.-Y. Le Traon, A. Pascual, A. Ponte, B. Qiu, Rascle N, C. Ubelmann, J. Wang, and E.D. Zaron.  Global Observations of Fine-Scale Ocean Surface Topography with the Surface Water and Ocean topography (SWOT) Mission.  2019\.
  * [33] D.R. Munday, H.L. Johnson, and D.P. Marshall.  Impacts and effects of mesoscale ocean eddies on ocean carbon storage and atmospheric pCO2.  Global Biogeochemical Cycles, 28(8):877–896, 2014.
  * [34] T. Palmer.  Build high-resolution global climate models.  Nature, 515:338–339, 2014.
  * [35] A. Ramadhan, G. Wagner, C. Hill, J.-M. Campin, V. Churavy, T. Besard, A. Souza, A. Edelman, R. Ferrari, and J. Marshall.  Oceananigans.jl: Fast and friendly geophysical fluid dynamics on GPUs.  Journal of Open Source Software, 5(53), 2020.
  * [36] T. Ringler, M. Petersen, R.L. Higdon, D. Jacobsen, P.W. Jones, and M. Maltrud.  A multi-resolution approach to global ocean modeling.  Ocean Modelling, 69:211–232, 2013.
  * [37] M.J. Roberts, A. Baker, E.W. Blockley, D. Calvert, A. Coward, H.T. Hewitt, L.C. Jackson, T. Kuhlbrodt, P. Mathiot, C.D. Roberts, R. Schiemann, J. Seddon, B. Vannière, and P.L. Vidale.  Description of the resolution hierarchy of the global coupled hadgem3-gc3.1 model as used in cmip6 highresmip experiments.  Geosci. Model Dev., 12(12):4999–5028, 2019.
  * [38] V. Churavy S. Byrne, L.C. Wilcox.  MPI.jl: Julia bindings for the message passing interface.  Proceedings of the JuliaCon Conferences, 2021.
  * [39] O.A. Saenko, D. Yang, and J.M. Gregory.  Impact of mesoscale eddy transfer on heat uptake in an eddy-parameterizing ocean model.  Journal of Climate, 31(20):8589 – 8606, 2018.
  * [40] T. Schneider, S. Behera, G. Boccaletti, C. Deser, K. Emanuel, R. Ferrari, L.R. Leung, N. Lin, T. Müller, A. Navarra, O. Ndiaye, A. Stuart, J. Tribbia, and T. Yamagata.  Harnessing AI, Data, and Computing to Advance Climate Modeling and Prediction.  Nature Climate Change, submitted.
  * [41] T. Schneider, S. Lan, A. Stuart, and J. Teixeira.  Earth System Modeling 2.0: A Blueprint for Models That Learn From Observations and Targeted High-Resolution Simulations.  Geophysical Research Letters, 2017.
  * [42] C. Shu.  Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws.  Icase report no. 97-65, Institute for Computer Applications in Science and Engineering, NASA Langley Research Center, 1997.
  * [43] S. Silvestri, G.L. Wagner, A. Souza, N. Constantinou, C. Hill, J.-M. Campin, and R. Ferrari.  A WENO-based Vector Invariant advection scheme for Implicit Ocean LES.  In preparation.
  * [44] D. Struck.  Meet the team shaking up climate models.  The Christian Science Monitor, 2021.
  * [45] Z. Su, J. Wang, P. Klein, A. Thompson, and D. Menemenlis.  Ocean submesoscales as a key component of the global heat budget.  Nat. Commun., 9:775, 2018.
  * [46] J.R. Taylor and A.F. Thompson.  Submesoscale Dynamics in the Upper Ocean.  Annual Review of Fluid Mechanics, 55(1):103–127, 2023.
  * [47] MPI Team.  MPI: A Message-Passing interface standard.  Technical report, USA, 1994.
  * [48] P. Wang, J. Jiang, P. Lin, M. Ding, J. Wei, F. Zhang, L. Zhao, Y. Li, Z. Yu, W. Zheng, Y. Yu, X. Chi, and H. Liu.  The GPU version of LASG/IAP climate system ocean model version 3 (LICOM3) under the heterogeneous-compute interface for portability (HIP) framework and its large-scale application.  Geosci. Model Dev., 14(5):2781–2799, 2021.
  * [49] N.P. Wedi, I. Polichtchouk, P. Dueben, V.G. Anantharaj, P. Bauer, S. Boussetta, P. Browne, W. Deconinck, W. Gaudin, I. Hadade, et al.  A baseline for global weather and climate simulations at 1 km resolution.  Journal of Advances in Modeling Earth Systems, 12(11):e2020MS002192, 2020.
  * [50] B. Xiao, F. Qiao, Q. Shu, X. Yin, G. Wang, and S. Wang.  The development and validation of a global 1/32° surface-wave-tide-circulation coupled ocean model: FIO-COM32.  Geosci. Model Dev., 2022.
  * [51] S. Zhang, H. Fu, L. Wu, Y. Li, H. Wang, Y. Zeng, X. Duan, W. Wan, L. Wang, Y. Zhuang, et al.  Optimizing high-resolution community earth system model on a heterogeneous many-core supercomputing platform.  Geoscientific Model Development, 13(10):4809–4829, 2020.


\regret^T(\mu^\star)&\leq H\sqrt{\iota'\logtwo\AX T}+2H\sqrt{2\iota'\logtwo \AX T}+2H\sqrt{\iota'\logtwo \AX T}\\
&\leq 6H\sqrt{\iota'\logtwo \AX T}
\end{align*}

where we used $2\sqrt{2}\leq 3$.


# Analyzing and comparing door-to-door travel times for air transportation
using aggregated Uber data

Philippe Monmousseau, Aude Marzuoli, Eric Feron and Daniel Delahaye

###### Abstract

Improving the passenger air travel experience is one of the explicit goals set
by the Next Generation Air Transportation System in the United States and by
the Advisory Council for Aeronautics Research in Europe FlightPath 2050. Both
suggest door-to-door travel times as a potential metric for these objectives.
In this paper, we propose a data-driven model to estimate door-to-door travel
times and compare the reach and performance of different access modes to a
city, as well as conduct segment analysis of full door-to-door trips. This
model can also be used to compare cities with respect to the integration of
their airport within their road structure. We showcase multiple applications
of this full door-to-door travel time model to demonstrate how the model can
be used to locate where progress can be made.

##  I Introduction

Both in Europe and in the United States, national or supra-national agencies
promote the need for seamless door-to-door travel and data sharing. They were
deemed as needed by the European Commission’s 2011 White Paper [1] and were
reconfirmed by the Federal Aviation Administration (FAA) in 2017 [2]. Data
sharing was already a main focus in the early 2000s; in response, Europe
created and adopted SWIM - System Wide Information Management [3] \- and the
FAA followed suit. The Next Generation Air Transportation System (NextGen) [4]
in the United States and the Advisory Council for Aeronautics Research in
Europe (ACARE) Flightpath 2050 [1] both aim to have a more passenger-centric
approach. To this end, ACARE Flightpath 2050 sets some ambitious goals, which
are not all measurable yet due to lack of available data. For example, it aims
at having 90% of travelers within Europe being able to complete their door-to-
door journey within 4 hours. In the US, the Joint Planning and Development
Office has proposed and tested metrics regarding NextGen’s goals [5], but the
passenger-centric metrics, especially regarding door-to-door travel times, are
still missing.

Cook et al. [6] first explored the shift from flight-centric to passenger-
centric metrics in the project POEM - Passenger Oriented Enhanced Metrics -
where they propose propagation-centric and passenger-oriented performance
metrics and compare them with existing flight-centric metrics. Later, Laplace
et al. [7] introduce the concept of Multimodal, Efficient Transportation in
Airports and Collaborative Decision Making (META-CDM); they propose to link
both airside CDM and landside CDM, thus taking passenger perspective into
account. In this perspective, Kim et al. [8] propose an airport gate
scheduling model for improved efficiency with a balance between aircraft,
operator and passenger objectives. Dray et al. [9] illustrate the importance
of multimodality by considering ground transportation as well during major
disturbances of the air transportation system in order to offer better
solutions to passengers.

The estimation of door-to-door travel time for multi-modal trips has been
previously studied, but for trips contained within the same metropolitan area.
Peer et al. [10] focus on commutes within a Dutch city by studying door-to-
door travel times and schedule delays for daily commuters, and show that, for
the estimation of the overall travel time, it is important to consider the
correlation of travel times across different road links. Salonen and Toivonen
[11] investigate the need for comparable models and measures for trips by car
or public transport with focus on the city of Helsinki. Their multi-modal
approach takes into account the walking and waiting times necessary to reach a
station or a parking place. Duran-Hormazabal and Tirachini [12] analyze travel
time variability for multi-modal trips within the city of Santiago, Chile,
using human surveyors and GPS data to estimate the time spent in the different
transportation modes, namely walking, driving a car, riding a bus and taking
the subway. Pels et al. [13] analyze the relative importance of access time to
airports in the passengers’ choice within the San Francisco Bay Area based on
a passenger survey, offering perspective from air transportation. These works
emphasize the importance of considering all relevant modes when estimating
door-to-door travel times, but are limited in scope with respect to the size
of the area considered and the amount of data available.

Thanks to the increasing use of mobile phones as data sources, larger scale
studies with a focus on air transportation have been possible. In the United
States, Marzuoli et al. [14] implement and validate a method to detect
domestic air passengers using mobile phone data available on a nationwide
scale. Though the main focus of this work is the passenger behavior at
airports, the granularity of the data facilitates analysis of each phase
within the full door-to-door trip. Marzuoli et al. [15] then combine mobile
phone data with social media data to analyze passenger experiences in airports
during a major disruptions. In Europe, within the BigData4ATM
project111www.bigdata4atm.eu, Garcia-Albertos et al. [16] also present a
methodology for measuring door-to-door travel times using mobile phone data,
illustrated through trips between Madrid and Barcelona. Mobile phone data are,
however, proprietary data, which are difficult to access for research.

Grimme and Martens [17] propose a model analyzing the feasibility of the
4-hour goal proposed by FlightPath 2050 based on airport-to-airport flight
times and a simplified model of access to and from airports. Sun et al. [18]
use open source maps and datasets to calculate door-to-door minimum travel
time estimations in order to study the possible competitiveness of air taxis.

In the upcoming sections, the model and analysis presented are also based on
already available online data but with a post operation approach. The aim of
this model is to create a method to measure the average door-to-door travel
times once trips are completed to analyze and compare available modes of
transportation. We have applied the first version of this method to two intra-
European multi-modal trips, thus comparing air transportation and rail
transportation [19]. We then used an improved version of the same method,
leveraging four different data sources (ride-sharing, flight, phone, and
census data) and adapted to the conditions in the United States, to compare
trips using direct flights between five US cities, three of them on the West
Coast and the other two on the East Coast [20].

In this paper, we offer a data-driven model for the computation of door-to-
door travel times that harnesses recently available data along with public
data. The data-driven methods developed can be applied for most multi-modal
trips between two cities where relevant data are available and are not limited
to the air transportation system. The range of new analyses available using
this model is illustrated with multiple modal analysis of an intra-European
trip, a per-leg analysis of multiple intra-USA trips and an analysis of the
impact of severe weather disruptions. These analyses have direct applications
for passengers, urban planners and decision makers and highlight the
difference between taking a flight-centric approach to the air transportation
system and taking a passenger-centric approach.

Section II of this paper presents the data-driven, full door-to-door travel
time model; Section III showcases a first set of analyses and applications
facilitated by this model for trips between Amsterdam and Paris; and, Section
IV focuses on a set of analyses for trips within the United States where more
data are available; Finally, Section V concludes this paper and proposes
future research directions.

##  II The full door-to-door data-driven model

Similarly to [16] and [18], we can deconstruct the travel time $T$ for trips
with direct flights or direct train links into five different trip phases,
represented in Figure 1 and summarized in equation (1),

$T=t_{\text{to}}+t_{\text{dep}}+t_{\text{in}}+t_{\text{arr}}+t_{\text{from}}\leavevmode\nobreak\
,$ (1)

where

  * •

$t_{\text{to}}$ is the time spent traveling from the start of the journey to
the departure station (e.g. train station or airport),

  * •

$t_{\text{dep}}$ is the time spent waiting and going through security
processes (if any) at the departure station,

  * •

$t_{\text{in}}$ is the time actually spent in flight or on rails,

  * •

$t_{\text{arr}}$ is the time spent at the arrival station (e.g. going through
security processes),

  * •

$t_{\text{from}}$ is the time spent traveling from the arrival station to the
final destination.

DeparturestationArrivalstation$t_{\text{to}}$$t_{\text{dep}}$$t_{\text{in}}$$t_{\text{arr}}$$t_{\text{from}}$
Figure 1: Model of the full door-to-door travel time.

The full model for door-to-door travel time proposed in this paper is
established by data-driven methods used to calculate the values of the
different times contained in equation (1). These data-driven methods are
described in Sections II-A through II-D.

This study focuses on air and rail transportation as main transportation
modes, which give the value of $t_{\text{in}}$, though the process can also be
applied to inter-city bus trips. Furthermore, it is assumed that passengers
travel by road when arriving or leaving the main station (airport or train
station) for the calculation of $t_{\text{to}}$ and $t_{\text{from}}$. In
response to data availability, the case studies only consider six major US
cities (Atlanta, Boston, Los Angeles, Seattle, San Francisco and Washington
D.C.) and two European capitals (Amsterdam and Paris).

###  II-A Travel time from the origin location to the departure station and
from the arrival station to the final destination

We can estimate the road transit times from origin location to departure
station ($t_{\text{to}}$) and from arrival station to final destination
($t_{\text{from}}$) by using aggregated and publicly available data from taxi
or ride-sharing services. Uber [21] is a ride-sharing service launched in 2010
and located in major urban areas on six continents; it has recently released
anonymized and aggregated travel time data for certain of the urban areas
where it operates. The available data consist of the average, minimum and
maximum travel times between different zones (e.g. census tracts in the case
of US cities) within serviced area from all Uber rides aggregated over five
different periods for each considered day. The five considered periods, used
throughout this study, are defined as follows:

  * •

Early Morning: from midnight to 7am

  * •

AM: from 7am to 10am

  * •

Midday: from 10am to 4pm

  * •

PM: from 4pm to 7pm

  * •

Late Evening: from 7pm to midnight

There are days when the travel times between some zones are only aggregated at
a daily level. Travel times are associated with their mean starting door time,
i.e. the mean of all the time stamps from the trip contained in the zone of
departure.

Since Uber was initially introduced in the US, the impact of Uber in US urban
transit has already been the focus of several studies prior to this data
release. Li et al. [22] concludes that, at an aggregated level, Uber tends to
decrease congestion in the US urban areas where it was introduced. Later,
Erhardt et al. [23] build a model showing that ride sharing companies do
increase congestion using the example of San Francisco. Hall et al. [24] focus
on whether Uber complemented or substituted public transit by studying the use
of public transit system before and after Uber’s entry date in different US
cities. Wang and Mu [25] study Uber’s accessibility in Atlanta, GA (US) by
using the average wait time for a ride as a proxy and conclude that the Uber
use is not associated to a specific social category. Following the release of
Uber data, Pearson et al. [26] propose a traffic flow model based on this
aggregated Uber data and use it to analyze traffic patterns for seven cities
world-wide. Assuming Uber rides as part of the road traffic flow, this study
considers that Uber’s travel times are an acceptable proxy of the actual
travel times by road. In cities where busses don’t have specific road lanes,
these travel times are a valid proxy for both car and bus trips. This paper
limits its scope to road access to and egress from considered stations. The
analysis of subway alternatives is not considered in this paper.

Each US city is divided into their census tracts; Paris into the IRIS zones
used by INSEE [27] for census, and Amsterdam into its official districts
called _wijk_.

###  II-B Dwell time at stations

The dwell time at a station, either $t_{\text{dep}}$ or $t_{\text{arr}}$, is
defined as the time spent at the station, whether going through security
processes, walking through the station, or waiting. The time spent at each
station depends on the mode considered, the specific trip, and whether the
passenger is departing or arriving. The dwell time at departure can be split
into two components,

$t_{\text{dep}}=t_{\text{sec}}+t_{\text{wait}}.\leavevmode\nobreak\ ,$ (2)

a processing time, $t_{\text{sec}}$, necessary to get through security (if
any) and through the station to the desired gate or track, and an extra wait
time, $t_{\text{wait}}$, due to unanticipated delays.

Processing times at US airports are based on the average wait times at
airports extracted from the study of Marzuoli et al. [14]. The six US airports
under study in this paper are: Hartsfield-Jackson Atlanta International
Airport (ATL), Boston’s Logan International Airport (BOS) and Ronald Reagan
Washington National Airport (DCA) for the East Coast, Los Angeles
International Airport (LAX), Seattle-Tacoma International Airport (SEA) and
San Francisco International Airport (SFO) for the West Coast. Processing times
at European airports are assumed invariant between airports and determined
using most airline recommendations. The three European airports under study
are: Paris Charles de Gaulle Airport (CDG), Paris Orly Airport (ORY) and
Amsterdam Airport Schiphol (AMS).

The average dwell times at these airports are summarized in Table II for US
airports and in Table II for European airports.

Table I: Average dwell time spent at US airports in minutes.

| ATL | BOS | DCA | LAX | SEA | SFO
---|---|---|---|---|---|---
Time at departure | 110 | 105 | 100 | 125 | 105 | 105
Time at arrival | 60 | 40 | 35 | 65 | 50 | 45

Table II: Average dwell time spent at European airports in minutes.

| AMS | CDG | ORY
---|---|---|---
Time at departure | 90 | 90 | 90
Time at arrival | 45 | 45 | 45

With regard to processing times at train stations, based on the recommendation
of the train station websites, the departure dwell time is set at 15 minutes
and the arrival dwell time is set at 10 minutes for all train stations. We can
improve these estimates by gathering data from GPS or mobile phone sources as
well as WiFi beacons within airports and train stations, and by using a method
similar to Nikoue et al. [28].

We can calculate the extra wait times when the scheduled and actual departure
or arrival times are available. For US airports, these wait times are
calculated only for departure using the publicly available data from the
Bureau of Transportation Statistics (BTS) [29]. They were obtained by
subtracting the scheduled departure time from the actual flight departure
time.

###  II-C Time in flight or on rail

####  II-C1 US flights

The actual flight time was calculated based on the data from BTS using the
actual departure/arrival times of all direct flights between each city pairs
from January 1${}^{\text{st}}$ 2018 to March 31${}^{\text{st}}$ 2018\.
Cancelled flights are not considered in this study and were discarded.

####  II-C2 European trips

In Europe, we assume that flights and trains are on time and follow a weekly
schedule, due to a lack of publicly centralized flight schedule data. The
weekly schedules are extracted from actual train and flight schedules gathered
over a period of several months and are assumed applicable over the full
period under study.

###  II-D Full door-to-door travel time

Our model assumes that travelers plan their departure time to arrive at the
departure station exactly $t_{\text{sec}}$ minutes (eq. (2)) before the
scheduled departure time of their flight or train. We use this assumption to
determine the value of $t_{\text{to}}$ since it defines the period of the day
to consider when extracting the Uber average time from the origin location to
the departure station. We extract the value of $t_{\text{from}}$ by using the
actual arrival time of the flight or train. When only daily aggregated times
are available in the Uber data, these times are used for each period of the
day in proxy.

##  III Flights versus trains: a comparison of different access modes to Paris

Let us consider a traveler leaving from Amsterdam city center to reach the
Paris area. We have chosen the city center of Amsterdam as it covers both
tourists and business travelers, but the proposed door-to-door travel time
model and subsequent analysis can be applied from any zone. Three possible
means of transportation are under study for this trip: plane from Amsterdam
Airport Schiphol to Paris Charles De Gaulle airport (CDG) or via Paris Orly
(ORY), or train from Amsterdam Centraal via Paris Gare du Nord (GDN).

###  III-A Flight and train schedules

As in Section II-C2, the flight and train schedules used for this study are
weekly schedules. The weekly flight schedules between Amsterdam Airport
Schiphol (AMS) and CDG or ORY were extracted from the actual flight schedules
from December 2019 to January 2020, and it was assumed that the obtained
weekly schedules would run from January 1${}^{\text{st}}$ 2018 to September
30${}^{\text{th}}$ 2019\. These weekly schedules are summarized in Table III-A
for flights between AMS and ORY, and in Table III-A for flights between AMS
and CDG.

Table III: Extracted weekly schedule from Amsterdam to Paris via ORY.

[table head=Mo Tu We Th Fr Sa Su Ams. Paris
]tables/ams2ory.csv  Table IV: Extracted weekly schedule from Amsterdam to
Paris via CDG. [table head=Mo Tu We Th Fr Sa Su Ams. Paris
]tables/ams2cdg.csv

The weekly train schedule between Amsterdam Centraal station and Paris Gare du
Nord (GDN) is similarly extracted from the actual train schedule of the year
2019 and applied to the year 2018. It is summarized in Table III-A. Night
trains are not considered for this study.

Table V: Extracted weekly schedule from Amsterdam to Paris via GDN.

[table head=Mo Tu We Th Fr Sa Su Ams. Paris
]tables/amc2gdn.csv

These schedules already highlight the major differences between the three
considered modes. The flight schedule through ORY contains the fewest
possibilities, limited to two flights daily, whereas the other two models
offer hourly scheduled transport. Another notable difference can be seen with
respect to the station-to-station travel times: flights between Amsterdam and
Paris (both CDG and ORY) take 1h20 ($\pm$ 5 minutes) while train rides between
Amsterdam Centraal and Paris GDN take 3h20 ($\pm$ 3 minutes).

###  III-B Average total travel time mode comparison

The proposed data-driven door-to-door model can be used to evaluate and
compare the range of each considered mode, which helps to understand better
the urban structure and behavior from a transportation point of view.

For each trip $\tau$ (flight or rail) over the period $\mathcal{D}$ from
January 1${}^{\text{st}}$ 2018 to September 30${}^{\text{th}}$ 2019, the
associated average full door-to-door travel time $\bar{T}(\tau)$,

$\bar{T}(\tau)=\bar{t}_{\text{to}}(\tau)+t_{\text{dep}}(\tau)+t_{\text{in}}(\tau)+t_{\text{arr}}(\tau)+\bar{t}_{\text{from}}(\tau)\leavevmode\nobreak\
,$ (3)

where $\bar{t}_{\text{to}}(\tau)$ is the average ride time between Amsterdam
city center and the departure station (AMS or Amsterdam Centraal) for the trip
$\tau$, $\bar{t}_{\text{from}}(\tau)$ is the average ride time from the
arriving station (CDG, ORY or GDN) to the final arrival zone for the trip
$\tau$.

The same daily periods as those used in the Uber data (see Section II-A) are
considered here to categorize the trips into five groups depending on the time
of arrival at the final destination. For each day $d$ and each period $p$, the
mean per arrival zone $z$ of the average door-to-door travel times is
calculated for each mode $m$,

$E^{d,p}_{z}(m)=\frac{1}{|\mathcal{T}^{d,p}_{z}|}\sum_{\tau\in\mathcal{T}^{d,p}_{z}}\bar{T}(\tau)\leavevmode\nobreak\
,$ (4)

where $\mathcal{T}^{d,p}_{z}$ is the set of flight and rail trips that end at
zone $z$ on day $d$ and period $p$.

For each day $d$, each period $p$ and each arrival zone $z$, the mode with the
shortest mean travel time $E^{d,p}_{z}(m)$ is kept. The number of times
$N^{p}_{z}(m)$ a mode $m$ has had the shortest mean travel time is counted for
each zone $z$ and for each period $p$ over the twenty-one month period
$\mathcal{D}$,

$N^{p}_{z}(m)=|\\{d\in\mathcal{D}\leavevmode\nobreak\ |\leavevmode\nobreak\
m=\text{arg}\min_{m_{1}}E^{d,p}_{z}(m_{1})\\}|\leavevmode\nobreak\ .$ (5)

This distribution of modes over the different zones can help travelers choose
the mode of transport that is best suited depending on the desired arrival
zone and on the desired time of arrival. It can also help urban planners to
better understand the road network linking the different stations to the city.

Figure 2 shows the fastest mode to reach the different zones in the Paris
dataset for the five different periods of the days used by the Uber dataset.
For each zone $z$ and each period $p$, the fastest mode associated is the mode
$m$ having the highest $N^{p}_{z}(m)$, i.e. the highest number of days with
the lowest average total travel time over the considered date range. The zones
best reached through CDG are indicated in blue, ORY in red and GDN in green.

(a) Morning (AM)

(b) Midday

Figure 2: Comparison of the average total travel times to the Paris area
between the three considered arrival stations (CDG: blue, ORY: red, GDN:
green) for a trip starting from Amsterdam city center for different trip
termination periods.

We can draw several conclusions from these maps. The absence of zones best
reached by plane via ORY (in red) is particularly noticeable in the morning
period: An important area of South-West Paris is not reached by Uber rides
neither from GDN nor from CDG. These maps would advocate for an increase in
frequency for the AMS-ORY flights from a traveler’s perspective.

From a structural perspective, the interstate linking Paris to CDG is visible
on all three maps since it enables travelers through GDN to reach zones close
to CDG faster than if they flew to CDG directly. The perimeter highway
circling Paris is also a major aid to GDN and is visible on the maps where
there is an important competition between GDN and CDG. The section of the
perimeter highway farthest from GDN (i.e. in the south-west of Paris) is,
however, overtaken by either airport depending on the period of the day. The
rest of GDN influence zone is fairly invariant from a period to another.

Using a similar map representation for trips ending in the afternoon but not
shown here for space considerations, ORY’s range is limited during the
afternoon, with CDG taking over some zones close to ORY. This is essentially
due to the limited number of flights landing in the afternoon (one per week,
every Friday) compared to the daily arrival of CDG flights.

###  III-C Average total travel time distribution analysis

Once the fastest mode to reach each zone is determined, it is possible to
analyze the fastest full door-to-door travel time for each zone. This approach
gives an overview of the level of integration of airports, train stations, and
road structure and can indicate zones that are less reachable than others and
would thus require more attention from urban planners.

The fastest full door-to-door travel time associated with a zone $z$ at a
period $p$ is calculated as the average fastest travel time to reach zone $z$
at period $p$ across all modes and over the period $\mathcal{D}$,

$\bar{E}_{z}^{p}=\frac{1}{|\mathcal{D}|}\sum\min_{m}E_{z}^{d,p}(m)\leavevmode\nobreak\
.$ (6)

Figure 3 displays the fastest full door-to-door travel times
$\\{\bar{E}_{z}^{p}\\}_{p,z}$ to reach the different zones in the Paris
dataset for trips finishing in the morning or at midday. The color scale that
indicates the fastest full door-to-door travel times is identical in all
subfigures. The contour of each zone indicates the fastest mode to reach it
according to the results from Section III-B using the same color code as
Figure 2, i.e. the zones reached faster through CDG are surrounded in blue,
ORY in red and GDN in green.

(a) Morning (AM)

(b) Midday

Figure 3: Comparison of the fastest full door-to-door travel times to the
Paris area between the three considered arrival stations (CDG: blue, ORY: red,
GDN: green) for a trip starting from Amsterdam city center for different trip
termination periods. The contour color of each zone indicates the fastest mode
to reach it.

For a better comparison, the distribution of the number of zones per period
reached within four time intervals (less than 4 hours, between 4 hours and 4
hours and 30 minutes, between 4 hours and 30 minutes and 5 hours, and more
than 5 hours) is presented in Table III-C.

Table VI: Number of zones per mode and period of the day grouped by full door-
to-door travel time intervals. The original dataset is the same as that used
to generate Figure 3.

Mode | Time interval | Early | AM | Midday | PM | Late
---|---|---|---|---|---|---
CDG |  $t\leq$ 4h | 0 | 4 | 6 | 5 | 11
4h $<t\leq$ 4h30 | 22 | 1306 | 866 | 1189 | 845
4h30 $<t\leq$ 5h | 0 | 653 | 433 | 498 | 384
$t>$ 5h | 187 | 0 | 15 | 13 | 11
GDN |  $t\leq$ 4h | 398 | 247 | 247 | 290 | 314
4h $<t\leq$ 4h30 | 775 | 818 | 731 | 719 | 641
4h30 $<t\leq$ 5h | 0 | 14 | 8 | 6 | 8
$t>$ 5h | 0 | 0 | 0 | 0 | 0
ORY |  $t\leq$ 4h | 0 | 0 | 0 | 0 | 0
4h $<t\leq$ 4h30 | 0 | 0 | 906 | 563 | 997
4h30 $<t\leq$ 5h | 0 | 0 | 397 | 259 | 425
$t>$ 5h | 0 | 0 | 0 | 0 | 1

We see a dissymmetry between the north and the south of Paris by looking only
at the time color scale in Figure 3. The number of green zones, i.e. zones
reachable in less than 4 hours and 20 minutes, and the surface covered by
these green zones are more important in the northern part of Paris than in the
southern part of Paris, including in areas close to Paris Orly airport.
Combining this observation with the contour color of each zone suggests that
Paris Orly airport is not as well integrated in the Parisian road structure as
Paris Charles de Gaulle airport, which can make it less attractive for
travelers desiring to travel by air and thus less competitive.

We can complete a more quantitative analysis from Table III-C, with some of
the main findings listed here:

  * •

Only 10% of the arrival zones are reached in less than 4 hours from Amsterdam
city center.

  * •

Zones that can be reached in less than 4 hours from Amsterdam city center are
overwhelmingly reached by train through Paris Gare du Nord (98%).

  * •

A trip from Amsterdam city center to Paris going through ORY always takes more
than 4 hours.

  * •

78% of the arrival zones are reached in less than 4 hours and 30 minutes from
Amsterdam city center when combining all three possible modes.

Therefore, we can use these results to assess how well the 4-hour goal from
ACARE FlightPath 2050 is engaged.

###  III-D Reliability issues

The proposed full door-to-door travel time model assumes that passengers
choose their departure time in order to arrive exactly $t_{\text{sec}}$ before
the scheduled departure of their plane or train and that they also know how
long it takes to reach the departure station. However, in reality, there is an
uncertainty in the time the traveler will spend reaching the airport and in
airport processing times. This uncertainty often leads to an additional buffer
time implying an earlier departure time for the traveler. Using the presented
model with the available data, we can find the most reliable mode to use per
arrival zone. The most reliable mode for a given arrival zone is defined as
the mode with the lowest variability in travel time, i.e. the mode where the
difference between the maximum travel time and the minimum travel time to
reach that zone is the lowest. This comparison is useful for passengers or
trips that require an accurate arrival time rather than a minimum travel time.

Figure 4 shows the most reliable mode on average to reach the different zones
in the Paris dataset for trips finishing in the morning or at midday. As for
the previous analysis, the period was determined using the departure time of
the full door-to-door trip and uses the same color code, i.e. the zones
reached most reliably through CDG are indicated in blue, ORY in red and GDN in
green. For each zone and each period of the day, the most reliable mode
associated is the mode having the highest number of days with the lowest
average variability travel time over the considered date range.

(a) Morning (AM)

(b) Midday

Figure 4: Comparison of the average variability of travel times to the Paris
area between the three considered arrival stations (CDG: blue, ORY: red, GDN:
green) for a trip starting from Amsterdam city center for different trip
termination periods.

Though Figure 4 and Figure 2 are similar, there are some major differences
between average efficiency and average reliability. For example, though it is
on average faster to reach the zones close to the highway leading to CDG by
train, after 10:00 it is safer from a time variability perspective to reach
them via CDG. From a reliability perspective, CDG has claimed the quasi
totality of the zones surrounding it, except in the early morning where trips
through GDN are still better. When we compare all three modes using this
metric, it appears that GDN has the greatest decrease in competitiveness, with
its range smaller than when considering average travel times.

###  III-E Impact of faster processing times

The major difference between air and rail travel is the necessary processing
time both at departure and at arrival. In this particular study, with a flight
time of about 80 minutes, the current assumption of a departure processing
time of 90 minutes implies that travelers spend more time at their departure
airport than in flight, which greatly impacts the rapidity of air travel. With
the presented model, one can modify these assumed processing times in order to
study the impact of improving these times both from an airport perspective and
a passenger perspective. Let’s assume that the processing time at airports is
improved from 90 to 60 minutes at departure and from 45 to 30 minutes at
arrival. These modifications could be achieved in reality considering that
this is an intra-Schengen trip and that there isn’t any border control.

(a) Early morning with faster processing times

(b) Early morning with normal processing times

Figure 5: Comparison of the average total travel times to the Paris area
assuming faster airport processing times between the three considered arrival
stations (CDG: blue, ORY: red, GDN: green) for a trip starting from Amsterdam
city center for trips arriving at destination early in the morning. The
corresponding map with normal processing times is also reproduced.

Figures 5& 6 show which is the fastest mode on average to reach the different
zones in the Paris dataset for trips arriving at destination early in the
morning or at midday.

The first major difference with this processing time improvement can be seen
for trips arriving in the early morning (Figure 5): all zones previously
reached through CDG are no longer accessed at this period since they were
associated to the 21:45 flight of the previous day. This indicates that all
trips through CDG start and end on the same day, with no trips finishing after
midnight.

(a) Midday with faster processing times

(b) Midday with normal processing times

Figure 6: Comparison of the average total travel times to the Paris area
assuming faster airport processing times between the three considered arrival
stations (CDG: blue, ORY: red, GDN: green) for a trip starting from Amsterdam
city center for trips arriving at destination at midday. The corresponding
maps with normal processing times are also reproduced.

Looking at trips arriving at midday (Figure 6), trips through CDG are greatly
advantaged by this time improvement, with CDG taking over more than half of
GDN’s previous influence zone. This range increase from CDG can be explained
both by faster door-to-door travel times and by the increase of trips through
CDG arriving at midday (rather than in the afternoon). As a matter of fact,
besides in the early morning, GDN loses its competitiveness against both
airports, with its range greatly shrinking in size. The competition between
CDG and ORY remains unchanged, which is understandable since they both
received the same processing time improvement.

A quantitative analysis similar to the analysis presented in Table III-C
concludes that all trips are now conducted in less than five hours and that
99.8% of the zones reachable are reached in less than 4h30. ORY sees some
major improvements with 97.5% of the zones best reached through it reached in
less than four hours (compared to no trips in less than 4h in the initial
model), while increasing the number of zones it reaches the fastest.

Using a map representation similar to Section III-C, but not presented here
due to space considerations, it is possible to notice a 20-30 minutes shift in
the time distribution for every period except for early morning trips since
train processing times were unchanged. The upper bound travel time is also
unchanged for trips arriving in the morning, which would indicate that for
some zones, the processing time improvement resulted in no improvement or even
a worsening of the full trip travel time. Besides that exception, in this case
a 45 minutes improvement in airport processing time leads only to a maximum of
30 minutes of average total travel time improvement due to the influence of
train trips through GDN.

##  IV A multi-modal analysis of the US air transportation system

Additional insights are gained from this full door-to-door travel time model
thanks to the availability of complementary data. The United States is a
federal state the size of a continent, therefore various aggregated and
centralized datasets are more easily available to all. Several of these
datasets are used in this section to add applications to the presented full
door-to-door model. This US study is limited to the period from January 1st
2018 to March 31st 2018.

###  IV-A Flight schedule

As presented in the model definition in Section II-C1, both the scheduled
flight times and the actual flight schedules of most domestic flights can be
obtained via the Department of Transportation Bureau of Transportation
Statistics (BTS) [29]. This study considers only the six US airports presented
in Section II-B, three East-coast airports - Hartsfield-Jackson Atlanta
International Airport (ATL), Boston’s Logan International Airport (BOS) and
Ronald Reagan Washington National Airport (DCA) - and three West-coast
airports - Los Angeles International Airport (LAX), Seattle-Tacoma
International Airport (SEA) and San Francisco International Airport (SFO).

During this three-month period, 38,826 flights were considered, which
corresponds to 3,523 early flights, 8,170 morning flights, 13,451 midday
flights, 6.695 afternoon flights and 6,987 evening flights. The full door-to-
door travel times were then calculated for each scheduled flight from January
1${}^{\text{st}}$, 2018 to March 31${}^{\text{st}}$, 2018, using the model
presented in Section II.

###  IV-B Leg analysis

We can use the full door-to-door model to better understand the time spent in
each leg proportionally to the time spent on the overall trip. For each trip,
we calculate the percentage of time spent at each phase based on the full
door-to-door travel time. We then calculate the average percentage time spent
for each phase and for each city pair trip. Figure 7 shows the bar plot of
these average percentage times for the thirty considered city pairs. The city
pairs are sorted according to the percentage of time spent in the actual
flight phase.

Figure 7: Bar plot of the average proportion of the time spent within each
phase of the full door-to-door journey for all thirty considered trips.

With the proposed full door-to-door model, for all considered trips,
passengers spend on average more time at the departure airport than riding to
and from the airports. This figure also shows that, with this model, for some
short-haul flights, such as between SFO and LAX or between BOS and DCA,
passengers spend on average more time at the departure airport than in the
plane. Refining the full door-to-door model by considering tailored airport
processing times $t_{\text{sec}}$ at departure depending on the city pair and
not only on the departure airport could lead to a different conclusion.
However, this modification of the model would require more access to passenger
data.

###  IV-C Airport integration comparison

The proposed full door-to-door model allows us to compare each airport’s
integration within its metropolitan area. Each census tract is associated with
an internal point within its boundaries, and this internal point can be used
to automatically calculate the distance between airports and each census tract
of their metropolitan area. The internal points were defined using an online
database222www.usboundary.com based on the US government 2010 census. Figure 8
shows the scatter plot of the average daily ride time to each airport versus
the geodesic distance to the airport for the six considered airport. The
geodesic distance is the shortest distance between two points along the
surface of the Earth. Additionally, the plot also figures a linear regression
of these average time with respect to the distance to the airport. A steeper
slope for the linear regression indicates that it takes longer to reach the
airport from a given distance.

Figure 8: Scatter plot of the average ride time to the airport $t_{\text{to}}$
versus the distance to the airport from January 1${}^{\text{st}}$ 2018 to
March 31${}^{\text{st}}$ 2018\. Straight lines indicate the linear regression
fit for each city.

Figure 8 highlights the disparity between the range of each airport within
available data: DCA has a range limited to 20 km while SFO attracts Uber
riders from more than 120 km away. The other four airports have a similar
range. The difference in slope of their associated linear regression is
nonetheless useful to rank their integration within their region of
attraction. From this perspective, Seattle has the best integrated airport,
i.e. the smallest slope, followed by Atlanta, Boston and then Los Angeles.

###  IV-D Impact of severe weather analysis

Using the same door-to-door travel time visualization process and applying it
to different days can be a tool to better analyze the effects of severe
weather perturbations on the full door-to-door journey. As an example, the
winter storm previously studied in [15] is analyzed for trips between
Washington D.C. and Boston using Figure 9. This winter storm hit the East
Coast of the United States on January 4${}^{\text{th}}$ 2018, and led to the
closure of two airports in New York City, along with the cancellation of the
majority of flights flying to or from the North-Eastern US coast.

(a) Before landfall, on January 2${}^{\text{nd}}$, 2018

(b) After landfall, on January 5${}^{\text{th}}$, 2018

Figure 9: Average door-to-door travel times from Washington D.C. city hall to
Boston over a single day, before and after the Bomb Cyclone of January 2018.

Figure 9 shows the map of the average full door-to-door travel times to reach
the Boston area starting from Washington D.C. city hall on January
2${}^{\text{nd}}$ 2018, before the landfall of this winter storm, and on
January 5${}^{\text{th}}$ 2018, after the landfall of the winter storm.

The color scales representing the full door-to-door travel time are identical
from one map to another. A shift towards the red is visible from January
2${}^{\text{nd}}$ 2018 (Figure 8(a)) to January 5${}^{\text{th}}$ 2018 (Figure
8(b)), along with some census tracts disappearing from the considered range on
January 5${}^{\text{th}}$ 2018, due to lack of sufficient Uber ride data.
These two observations indicate that the full door-to-door travel times are
closer to the maximum average travel time than from the minimum travel time on
January 5${}^{\text{th}}$ 2018 compared to January 2${}^{\text{nd}}$ 2018, and
that some zones might have been sufficiently adversely impacted by the weather
to prevent rides from the airport to reach them.

###  IV-E On the importance of a passenger-centric approach to delays

A final application to the full door-to-door model presented in this paper
emphasizes the difference between flight delay and passenger delay. Since Uber
splits the day into five different periods, each with their traffic
idiosyncrasies with respect to peak times, we can calculate how much extra
travel time is required for a passenger when a flight does not arrive in the
scheduled period. For example, a flight expected to arrive in the early
morning that lands after 10:00 AM could result in the passenger getting
stranded in traffic when trying to leave the airport. Though airlines are not
responsible for road traffic, passengers can choose flights based on their
arrival time to avoid peak time traffic.

To calculate this extra travel at aggregated level, we have calculated the
difference of average travel time between the two periods concerned by flights
not arriving according to schedule for each arrival zone. These travel time
differences are then aggregated into one travel time difference per city pair
by weighing the travel time associated to each census tract with the
proportion of passengers initiating their trips from there, or finishing their
trip there. The number of passengers originating from or finishing within a
census tract is assumed to be proportional to the population density of the
considered census tract.

Another measure of sensitivity is to consider the maximum difference between
the maximum travel times of each zone between the two considered periods. This
second measure indicates the worst variation of the travel time upper bound,
i.e. the maximum difference a traveler can experience going from the airport
to their final destination zone.

Let us consider the flight UA460 from LAX to SFO scheduled to arrive on
Thursday February 15, 2018 at 18:02 local time and that landed with a minor
delay of 16 minutes. Due to the 45-minute processing time required to leave
the airport, this 16-minute delay shifts the departure period from the airport
from afternoon (PM) to late evening. The aggregated average extra travel time
is of 15 minutes and 40 seconds, i.e. a 16-minute flight delay triggered an
average 31-minute total delay for the passengers. Looking at the second
considered measure, the maximum travel time difference for this flight delay
is of 72 minutes, meaning that one passenger could experience a delay of 88
minutes resulting from this 16-minute flight delay. This first example
illustrates that passenger delays and aircraft delays are distinct.

Paradoxically, arriving earlier than scheduled for a flight does not
necessarily mean that the full door-to-door trip ends earlier. For example,
flight VX1929 from LAX to SFO scheduled to arrive on Thursday February 8, 2018
at 15:22 local time actually landed 25 minutes earlier. This implied that the
passengers were no longer leaving the airport in the afternoon (PM) period but
at midday. The aggregated average extra travel time is here of 15 minutes and
2 seconds, so on average travelers did arrive earlier than scheduled, but only
by about ten minutes and not the twenty-five minutes announced by the airline.
However, looking at the second measurement method again, the maximum ride time
difference is 66 minutes and 44 seconds, which means that a passenger could
end up arriving forty minutes later than if the flight had landed on time.

##  V Discussion & Conclusion

By leveraging Uber’s recently released data and combining it with several
other available data sources, this paper proposed a data-driven model for the
estimation of full door-to-door travel times for multi-modal trips both in
Europe and in the United States. Though the model is used for one city pair in
Europe and six different cities in the United States, it can be implemented
between any world city pair with sufficient available ride-sharing or taxi
data. The proposed model can be adapted depending on how much data about the
considered travel modes are available, since a weekly schedule containing no
information relative to delays can already lead to some meaningful insights
from both a passenger perspective and a city planner perspective.

Once aggregated at a city level, the presented door-to-door travel time model
can be used for a paired comparison of the different phases of the full door-
to-door journey between two cities. It also enables us to analyze the actual
time spent while travelling between two specific cities. The evaluation on a
national level of some of the passenger-centric objectives proposed within
NextGen in the United States and within ACARE FlightPath 2050 in Europe is
possible thanks to the proposed model, especially the objectives regarding how
well integrated airports are within their cities. The model can also provide
insight to how multi-modal trips are affected by severe weather disruptions,
indicating where improvements can be made. It also brings a valuable
measurement of the difference between flight delays and passenger delays,
emphasizing the need for passenger-centric metrics to evaluate the performance
of the air transportation system, which is not solely constituted of planes.

Future studies should consider integrating additional data, such as
alternative transportation modes, e.g. the subway, in the model once they are
available and when calculating the time needed to reach the departure station
or leave from the arrival station. Additionally, the knowledge of the actual
daily proportion of passengers travelling via the different approach modes
(road or rail) would lead to a better precision of the proposed full door-to-
door travel time model. Aggregated information from GPS or mobile phone
sources could possibly be used to determine this proportion without infringing
passenger privacy.

## Acknowledgments

The authors would like to thank Nikunj Oza from NASA-Ames, the BDAI team from
Verizon Media in Sunnyvale as well as the _Ecole Nationale de l’Aviation
Civile_ and King Abdullah University of Science and Technology for their
financial support. The authors are also grateful for the help of Marine
Lebeater for her feedback. The authors would also like to thank the SESAR
Joint Undertaking for its support of the project ER4-10 ”TRANSIT: Travel
Information Management for Seamless Intermodal Transport”. Data retrieved from
Uber Movement, (c) 2020 Uber Technologies, Inc., https://movement.uber.com.

## References

  * [1] M. Darecki, C. Edelstenne, E. Fernandez, P. Hartman, J.-P. Herteman, M. Kerkloch, I. King, P. Ky, M. Mathieu, G. Orsi, G. Schotman, C. Smith, and J.-D. Wörner, _Flightpath 2050: Europe’s Vision for Aviation ; Maintaining Global Leadership and Serving Society’s Needs ; Report of the High-Level Group on Aviation Research_ , E. Commission, Ed. Luxembourg: European Commission, 2011.
  * [2] _NextGen Priorities - Joint Implementation Plan Update Including the Northeast Corridor_. Federal Aviation Administration, Oct. 2017.
  * [3] J. Meserole and J. Moore, “What is System Wide Information Management (SWIM)?” in _2006 IEEE/AIAA 25TH Digital Avionics Systems Conference_. Portland, OR, USA: IEEE, Oct. 2006, pp. 1–8.
  * [4] NextGen Integration and Implementation Office, “NextGen Implementation Plan,” in _Federal Aviation Administration_ , 2009.
  * [5] Y. O. Gawdiak and T. Diana, “NextGen Metrics for the Joint Planning and Development Office,” vol. 11th AIAA Aviation Technology, Integration, and Operations (ATIO) Conference, including the AIAA Balloon Systems Conference and 19th AIAA Lighter-Than, 2011.
  * [6] A. Cook, G. Tanner, S. Cristóbal, and M. Zanin, “Passenger-Oriented Enhanced Metrics,” 2012.
  * [7] I. Laplace, A. Marzuoli, and E. Feron, “META-CDM: Multimodal, Efficient Transportation in Airports and Collaborative Decision Making,” in _Airports in Urban Networks 2014 (AUN2014)_ , Paris, 2014.
  * [8] S. H. Kim, A. Marzuoli, J.-P. Clarke, D. Delahaye, and E. Feron, “Airport Gate Scheduling for Passengers, Aircraft, and Operation,” in _Tenth USA/Europe Air Traffic Management Research and Development Seminar (ATM2013)_ , Chicago, Illinois, Jun. 2013.
  * [9] L. Dray, A. Marzuoli, and A. Evans, “Air Transportation and Multimodal, Collaborative Decision Making during Adverse Events,” in _Eleventh USA/Europe Air Traffic Management Research and Development Seminar (ATM2015)_ , Lisboa, Portugal, Jun. 2015.
  * [10] S. Peer, J. Knockaert, P. Koster, Y.-Y. Tseng, and E. T. Verhoef, “Door-to-door travel times in Revealed Preference departure time choice models: An approximation method using GPS data,” _Transportation Research Part B: Methodological_ , vol. 58, pp. 134–150, 2013.
  * [11] M. Salonen and T. Toivonen, “Modelling travel time in urban networks: Comparable measures for private car and public transport,” _Journal of Transport Geography_ , vol. 31, pp. 143–153, Jul. 2013.
  * [12] E. Durán-Hormazábal and A. Tirachini, “Estimation of travel time variability for cars, buses, metro and door-to-door public transport trips in Santiago, Chile,” _Research in Transportation Economics_ , vol. 59, pp. 26–39, Nov. 2016.
  * [13] E. Pels, P. Nijkamp, and P. Rietveld, “Access to and competition between airports: A case study for the San Francisco Bay area,” _Transportation Research Part A: Policy and Practice_ , vol. 37, no. 1, pp. 71–83, Jan. 2003.
  * [14] A. Marzuoli, E. Boidot, E. Feron, and A. Srivastava, “Implementing and Validating Air Passenger–Centric Metrics Using Mobile Phone Data,” _Journal of Aerospace Information Systems_ , vol. 16, no. 4, pp. 132–147, Apr. 2019.
  * [15] A. Marzuoli, P. Monmousseau, and E. Feron, “Passenger-centric metrics for Air Transportation leveraging mobile phone and Twitter data,” in _Data-Driven Intelligent Transportation Workshop - IEEE International Conference on Data Mining 2018_ , Singapore, Nov. 2018.
  * [16] P. García-Albertos, O. G. Cantú Ros, R. Herranz, and C. Ciruelos, “Understanding Door-to-Door Travel Times from Opportunistically Collected Mobile Phone Records,” in _SESAR Innovation Days 2017_ , 2017.
  * [17] W. Grimme and S. Maertens, “Flightpath 2050 revisited – An analysis of the 4-hour-goal using flight schedules and origin-destination passenger demand data,” _Transportation Research Procedia_ , vol. 43, pp. 147–155, 2019.
  * [18] X. Sun, S. Wandelt, and E. Stumpf, “Competitiveness of on-demand air taxis regarding door-to-door travel time: A race through Europe,” _Transportation Research Part E: Logistics and Transportation Review_ , vol. 119, pp. 1–18, Nov. 2018.
  * [19] P. Monmousseau, D. Delahaye, A. Marzuoli, and E. Feron, “Door-to-door travel time analysis from Paris to London and Amsterdam using Uber data,” in _Ninth SESAR Innovation Days_ , Athens, Greece, 2019.
  * [20] ——, “Door-to-door Air Travel Time Analysis in the United States using Uber Data,” in _1st International Conference on Artificial Intelligence and Data Analytics for Air Transportation (AIDA-AT 2020)_ , Singapore, Feb. 2020.
  * [21] Uber Technologies Inc., “Uber, About us,” 2020. [Online]. Available: http://www.uber.com/about
  * [22] Z. Li, Y. Hong, and Z. Zhang, “Do On-demand Ride-sharing Services Affect Traffic Congestion? Evidence from Uber Entry,” 2016.
  * [23] G. D. Erhardt, S. Roy, D. Cooper, B. Sana, M. Chen, and J. Castiglione, “Do transportation network companies decrease or increase congestion?” _Science Advances_ , vol. 5, no. 5, May 2019.
  * [24] J. D. Hall, C. Palsson, and J. Price, “Is Uber a substitute or complement for public transit?” _Journal of Urban Economics_ , vol. 108, pp. 36–50, Nov. 2018.
  * [25] M. Wang and L. Mu, “Spatial disparities of Uber accessibility: An exploratory analysis in Atlanta, USA,” _Computers, Environment and Urban Systems_ , vol. 67, pp. 169–175, Jan. 2018.
  * [26] M. Pearson, J. Sagastuy, and S. Samaniego, “Traffic Flow Analysis Using Uber Movement Data,” 2018.
  * [27] INSEE, “Institut national de la statistique et des études économiques,” 2020\. [Online]. Available: https://www.insee.fr
  * [28] H. Nikoue, A. Marzuoli, J.-P. Clarke, E. Feron, and J. Peters, “Passenger Flow Predictions at Sydney International Airport: A Data-Driven Queuing Approach,” _arXiv:1508.04839 [cs]_ , Aug. 2015.
  * [29] Bureau of Transportation Statistics, “Bureau of Transportation Statistics, About BTS,” 2018. [Online]. Available: http://www.rita.dot.gov/bts/about