Daniel Pereira
email: daniel.pereira.1@ens.etsmtl.ca

Guillaume Haiat

Julio Fernandes

Pierre Belanger

Effect of intracortical bone properties on the phase velocity and cut-off frequency of low-frequency guided wave modes (20-85 kHz)

The assessment of intracortical bone properties is of interest since early-stage osteoporosis is associated with resorption in the endosteal region. However, understanding the interaction between ultrasonic guided waves and the cortical bone structure remains challenging. The purpose of this work is to investigate the effect of intracortical bone properties on the ultrasonic response obtained at low-frequency (<100 kHz) using an axial transmission configuration. The semi-analytical finite element method was used to simulate the propagation of guided waves in a waveguide with realistic geometry and material properties. An array of 20 receivers was used to calculate the phase velocity and cut-off frequency of the excited modes using the 2D Fourier transform. The results show that the position of the emitter around the circumference of the bone is an important parameter to control since it can lead to variations of up to 10 dB in the amplitude of the transmitted modes. The cut-off frequency of the high order modes was, however, only slightly affected by the circunferential position of the emitter, and was sensitive mainly to the axial shear modulus. The phase velocity and cut-off frequency in the 20-85 kHz range are promising parameters for the assessment of intracortical properties.

INTRODUCTION

Osteoporosis is associated with a reduction of the bone mass and microarchitectural bone deterioration, leading to a reduction in bone quality [START_REF] Giangregorio | Fragility fractures and the osteoporosis care gap: An international phenomenon[END_REF][START_REF] Papaioannou | The osteoporosis care gap in canada[END_REF][START_REF] Ström | Osteoporosis: burden, health care provision and opportunities in the eu[END_REF] and to an increase in fracture risk (Consensus development conference: prophylaxis and treatment of osteoporosis, 1991). Endosteal resorption of cortical bone results in an increase in bone porosity at the inner part of the cortical shell and in a reduction of the cortical bone thickness [START_REF] Ritzel | The thickness of human vertebral cortical bone and its changes in aging and osteoporosis: A histomorphometric analysis of the complete spinal column from thirty-seven autopsy specimens[END_REF]. Dual-energy X-ray absorptiometry (DEXA), which is currently the gold standard for osteoporosis diagnostics [START_REF] Kanis | Assessment of fracture risk and its application to screening for postmenopausal osteoporosis: Synopsis of a who report[END_REF], nonetheless suffers from several limitations [START_REF] Gluer | A new quality of bone ultrasound research[END_REF][START_REF] Haba | Bone mineral densities and mechanical properties of retrieved femoral The Scientific World Journal[END_REF].

Quantitative ultrasound (QUS) was developed to assess bone quality and has the advantage of being non-radiative, non-invasive and relatively cheap [START_REF] Gluer | Quantitative ultrasound techniques for the assessment of osteoporosis: Expert agreement on current status[END_REF]. It, therefore, holds promise as a rapid screening method in a clinical setting. The mechanical nature of ultrasonic waves allows the technique to be used to retrieve the biomechanical properties of bone tissue [START_REF] Kaufman | Perspectives: Ultrasound assessment of bone[END_REF]. QUS techniques can be used to go beyond a simple estimation of the bone mineral density [START_REF] Nicholson | Ultrasound and the biomechanical competence of bone[END_REF], which can be achieved with DEXA and which is not sufficient to assess fracture risk (Office of the Surgeon General, 2004).

Bone QUS was initially developed in the context of trabecular bone characterization using transverse transmission devices [START_REF] Stein | Clinical assessment of the 1/3 radius using a new desktop ultrasonic bone densitometer[END_REF]. However, the investigation of cortical bone [START_REF] Rico | The therapy of osteoporosis and the importance of cortical bone[END_REF] has attracted significant interest since about 80% of the skeleton is made of cortical bone, which supports most of the body load, and is involved in osteoporotic fractures [START_REF] Seeley | Which fractures are associated with low appendicular bone mass in elderly women?[END_REF]. QUS using axial transmission techniques have mostly been used to study the cortical bone quality [START_REF] Haiat | Equivalent contributing depth investigated by a lateral wave with axial transmission in viscoelastic cortical bone[END_REF]. Cortical bones, such as the radius or the tibia, constitute suitable waveguides for the propagation of ultrasonic guided waves, as has was shown in various studies, including for instance (Gluer, 2008 , 1997; Haba et al., 2016). Ultrasonic guided waves have the advantage of being sensitive to both the mechanical and geometrical properties of cortical bone [START_REF] Foiret | Combined estimation of thickness and velocities using ultrasound guided waves: a pioneering study on in vitro cortical bone samples[END_REF][START_REF] Muller | Comparison of three ultrasonic axial transmission methods for bone assessment[END_REF][START_REF] Rozental | Premenopausal women with a distal radial fracture have deteriorated trabecular bone density and morphology compared with controls without a fracture[END_REF].

Most studies focusing on the assessment of cortical bone properties using axial transmission consider plate or cylinder waveguides (as an approximation of the actual bone geometry) to simplify the interpretation of ultrasonic responses obtained numerically or experimentally [START_REF] Kilappa | Tailoring the excitation of fundamental flexural guide waves in coated bone by phase-delayed array: Two-dimensional simulations[END_REF][START_REF] Le | Probing long bones with ultrasonic body waves[END_REF][START_REF] Minonzio | A free plate model can predict guided modes propagating in tubular bone-mimicking phantoms[END_REF][START_REF] Xu | Sparse svd method for high-resolution extraction of the dispersion curves of ultrasonic guided waves[END_REF]. Recently, more comprehensive numerical approaches were explored with the aim of improving the realism of the effect of cortical bone features on the axial transmission propagation expected in vivo. [START_REF] Bossy | Three-dimensional simulations of ultrasonic axial transmission velocity measurement on cortical bone models[END_REF] [START_REF] Bossy | Three-dimensional simulations of ultrasonic axial transmission velocity measurement on cortical bone models[END_REF] performed three-dimensional (3D) finite difference simulations on geometries derived from a human radius tomography image to evaluate the effect of bone curvature, anisotropy, and micro-porosity on the first arriving signal (FAS) velocity. However, the actual anatomical variations of bone biomechanical properties were not taken into account and more emphasis was brought to the analysis of numerical simulations performed with idealized 3D objects (tubes and semi-infinite hemicylinder geometries). [START_REF] Moilanen | Assessment of the cortical bone thickness using ultrasonic guided waves: Modelling and in vitro study[END_REF] [START_REF] Moilanen | Assessment of the cortical bone thickness using ultrasonic guided waves: Modelling and in vitro study[END_REF] developed a 2D numerical bone model to investigate the impact of realistic bone geometry specimens on the ultrasonic cortical thickness evaluation using plate assumptions. Haiat et al. (2009) [START_REF] Haïat | Influence of a gradient of material properties on ultrasonic wave propagation in cortical bone: Application to axial transmission[END_REF] [START_REF] Haiat | Equivalent contributing depth investigated by a lateral wave with axial transmission in viscoelastic cortical bone[END_REF] and [START_REF] Naili | Influence of viscoelastic and viscous absorption on ultrasonic wave propagation in cortical bone: Application to axial transmission[END_REF] [START_REF] Naili | Influence of viscoelastic and viscous absorption on ultrasonic wave propagation in cortical bone: Application to axial transmission[END_REF] assessed the effect of the heterogeneous nature of the cortical bone on the axial transmission response at 1 MHz in the framework of an anisotropic material by using 2D finite element (FE) simulation on a plate model. More recently, [START_REF] Chen | On ultrasound waves guided by bones with coupled soft tissues: A mechanism study and in vitro calibration[END_REF] [START_REF] Chen | On ultrasound waves guided by bones with coupled soft tissues: A mechanism study and in vitro calibration[END_REF] proposed a quantitative compensation for the effect of soft tissues based on an in vitro calibration, which facilitates the development of high-precision measurements of guided wave modes. Moreau et al. (2015) [START_REF] Moreau | Measuring the wavenumber of guided modes in waveguides with linearly varying thickness[END_REF] introduced a modified method to predict the dispersion curves of an isotropic plate waveguide with a linearly varying thickness along the propagation direction. The method allowed the detection of enhanced and more accurate wavenumbers in the context of the formulation of inverse problems. Recently, [START_REF] Tran | Sensitivity analysis of ultrasonic guided waves propagating in trilayered bone models: a numerical study[END_REF] [START_REF] Tran | Sensitivity analysis of ultrasonic guided waves propagating in trilayered bone models: a numerical study[END_REF] performed a numerical study to evaluate the effect of cortical thickness, stiffness coefficient, and thickness of the overlying soft tissues on the responses of fundamental ultrasonic guided waves in the frequency domain, using a homogeneous transversely isotropic tri-layered plate model. However, despite the improvements achieved so far concerning the modeling of cortical bone, a more detailed numerical study remains needed to explain the influence of the cross-sectional curvature and distribution of properties in the radial direction on the propagation of the guided waves at low-frequency.

Low-frequency (typically below 200 kHz) axial transmission research [START_REF] Egorov | Osteoporosis detection in postmenopausal women using axial transmission multi-frequency bone ultrasonometer: Clinical findings[END_REF][START_REF] Kilappa | Low-frequency axial ultrasound velocity correlates with bone mineral density and cortical thickness in the radius and tibia in pre-and postmenopausal women[END_REF][START_REF] Sarvazyan | Application of the dual-frequency ultrasonometer for osteoporosis detection[END_REF][START_REF] Tran | Sensitivity analysis of leaky lamb modes to the thickness and material properties of cor-tical bone with soft tissue: A semi-analytical finite element based simulation study[END_REF] demonstrated that ultrasonic guided waves are sensitive to changes in bone properties, such as the cortical thickness and the porosity. [START_REF] Muller | Comparison of three ultrasonic axial transmission methods for bone assessment[END_REF] [START_REF] Muller | Comparison of three ultrasonic axial transmission methods for bone assessment[END_REF] introduced low-frequency axial transmission as a promising method for assessing the cortical thickness. [START_REF] Tatarinov | Use of multiple acoustic wave modes for assessment of long bones: Model study[END_REF] [START_REF] Tatarinov | Use of multiple acoustic wave modes for assessment of long bones: Model study[END_REF] showed the potential of low-frequency (≈ 100 kHz) in assessing changes of deep underlying spongy layers in bovine tibia. More recently, Killappa et al.

(2015) [START_REF] Kilappa | Tailoring the excitation of fundamental flexural guide waves in coated bone by phase-delayed array: Two-dimensional simulations[END_REF] reported an improved performance of fundamental flexural modes used to assess cortical bone thickness using phase-delayed excitation at very lowfrequency (50 kHz). [START_REF] Kassou | Dispersion characteristics of the flexural wave assessed using low frequency (50-150khz) point-contact transducers: A feasibility study on bone-mimicking phantoms[END_REF] [START_REF] Kassou | Dispersion characteristics of the flexural wave assessed using low frequency (50-150khz) point-contact transducers: A feasibility study on bone-mimicking phantoms[END_REF] investigated the feasibility of dry point-contact transducers to infer the thickness of the in-vivo subjects in the 50-150 kHz frequency range. Operating at low-frequency offers the advantage of a lower attenuation as well as a reduced number of generated modes, which simplifies the signal analysis. In addition, low-frequency guided wave modes tend to achieve greater penetration depths due to the longer wavelengths as compared to high-frequency. As a consequence, the sensitivity to variations in intracortical bone properties could be improved when using low-frequencies, which is relevant to assess early stages of osteoporosis. Notwithstanding the remarkable level of realism that has been introduced in simulations over the past few years, the physical interaction between low-frequency guided waves and the cortical bone structure remains unexplored because the cortical bone has i) a complex cross-sectional geometry and ii) a heterogeneous distribution of material properties along the circumferential direction.

In a previous study by our group [START_REF] Pereira | Simulation of acoustic guided wave propagation in cortical bone using a semi-analytical finite element method[END_REF], the propagation of guided waves in an irregular, multi-layer and heterogeneous bone cross-section modeled with anisotropic and viscoelastic material properties was investigated. The effect of the intracortical bone properties was then evaluated using the first arriving signal (FAS) velocity obtained from only five receivers. However, due to the interaction between several modes around the FAS, the velocity associated with the FAS was shown to be a poor discriminator of the intracortical bone properties. The goal of this paper is therefore to determine the effect of intracortical bone properties on the phase velocity and cut-off frequency of low-frequency guided wave modes using an array of receivers, as well as a time-spatial frequency analysis technique instead of the FAS velocity. In order to do so, the excitation was generated by an emitter while the acquisition was performed in the time domain by an array of 20 receivers equally spaced along the bone surface. The bone ultrasonic responses were converted into the frequency domain using the two-dimensional Fourier transform (2D-FFT) in order to obtain the phase velocity and cut-off frequency of the propagating modes. Therefore, in this paper, the individual contribution of each viscoelastic coefficient and of the density was evaluated separately. The originality of this study when compared to our previous publication also lies in the application of the excitation source at different positions on the bone surface. This leads to a more comprehensive understanding of the physical interaction between the cortical bone curvature and the excitability of low-frequency guided wave modes.

The paper is organized as follows: section II.A introduces the semi-analytical finite element (SAFE) method used to simulate the signal of ultrasonic guided waves propagation in an arbitrary cross-section waveguide; section II.B describes the emitter/receivers configuration used to measure the ultrasonic responses on the bone surface in order to compute the frequency vs. wavenumber diagram. This section also describes the procedure employed to derive the phase velocity and cut-off frequency from the measured diagrams. In section II.C and II.D, the cross-sectional geometry of the cortical bone model and the methodology used to define the distribution of properties along the radial direction are presented respectively.

Section II.E focuses on the description of the experimental work performed on bone phan-toms that was done in order to verify the robustness and accuracy of the SAFE simulations.

In section III.A, the effect of the excitation position and the physiopathological condition on the phase velocity and cut-off frequency are presented. In addition, the individual contribution of each viscoelastic coefficient and density on the overall physiopathological conditions are presented in the section III.B. Conclusion including potential clinical implications are presented in section IV.

II. MATERIAL AND METHODS

A. Simulating time domain signals

The SAFE method has recently been used to simulate the modal properties of the ultrasonic guided wave in bone systems [START_REF] Nguyen | Semi-analytical solution of transient plane waves transmitted through a transversely isotropic poroelastic plate immersed in fluid[END_REF][START_REF] Nguyen | Computing dispersion curves of elastic/viscoelastic transversely-isotropic bone plates coupled with soft tissue and marrow using semi-analytical finite element (safe) method[END_REF][START_REF] Thakare | Semi-analytical finite-element modeling approach for guided wave assessment of mechanical degradation in bones[END_REF][START_REF] Tran | Sensitivity analysis of leaky lamb modes to the thickness and material properties of cor-tical bone with soft tissue: A semi-analytical finite element based simulation study[END_REF]. The computation of the time response from the modal properties of waveguides in an infinite arbitrary cross-section waveguide was previously described by [START_REF] Wilcox | Mode and transducer selection for long range lamb wave inspection[END_REF] and Loveday [START_REF] Loveday | Simulation of piezoelectric excitation of guided waves using waveguide finite elements[END_REF]. The resulting time domain signals consist in a linear superposition of the propagated response of each mode supported by the waveguide. A more detailed description of the implementation of the SAFE method can be found in the literature [START_REF] Fan | Torsional waves propagation along a waveguide of arbitrary cross section immersed in a perfect fluid[END_REF][START_REF] Predoi | Wave propagation along transversely periodic structures[END_REF], while the detailed procedure used to generate the time domain signal, can be found in a recently published study [START_REF] Pereira | Simulation of acoustic guided wave propagation in cortical bone using a semi-analytical finite element method[END_REF].

The dispersive response of a given mode excited from a point source in the out-of-plane direction (radial direction r in Fig. 1) can be calculated in the time domain using the following expression [START_REF] Wilcox | Mode and transducer selection for long range lamb wave inspection[END_REF]:

u (t) = 1 2π +∞ -∞ F (ω)E (ω) H (1) 0 (k (ω) z 0 ) e -iωt dω (1) i = √ -1 (2)
where u(t) is the surface out-of-plane displacement for a given mode as a function of time t, for an arbitrary propagation distance z 0 . The term H

(1) 0 is the zeroth-order Hankel function of the first kind and F (ω) is the frequency spectrum of a force input signal u input (t). The term k(ω) is the mode wavenumber, defined as a complex number by:

k(ω) = k(ω) real + ik(ω) imag (3)
where the real part denotes the propagating term and the imaginary part denotes the attenuation associated with the wavenumber. The term E(ω) is the frequency-dependent out-of-plane excitability of a given mode, and can be obtained with [START_REF] Wilcox | The excitation and detection of lamb waves with planar coil electromagnetic acoustic transducers[END_REF]:

E(ω) = ik(ω)ω 8 u out (ω) 2 P z (ω) (4)
where P z (ω) is the total power flow in the z-direction associated with the mode-shape as a function of the angular frequency, and u out (ω) is the out-of-plane displacement on the surface of the mode-shape (at the position of excitation) as a function of the angular frequency. The dispersion curves of the waveguide (k, P z and u out ) were obtained by solving the SAFE equations using the partial differential equation package in the Comsol Multiphysics and Matlab LiveLink environment. The complete time domain signal u sum (t) can then be obtained with the linear superposition of n modes supported by the waveguide by:

u sum (t) = n m=1 u m (t) (5)
where u m (t) is the surface out-of-plane displacement associated with the m-th mode and n is the total number of modes supported by the waveguide. The advantage of the SAFE method as compared to other simulation method such as 3D FE modeling is that with SAFE, only the cross-section is modeled and meshed, resulting in computationally efficient simulations, which may represent a significant gain of computation time. Furthermore, the dispersion curves are only calculated once using the SAFE method and then a number of time domain signals can be generated within minutes at minimal computational cost. The equivalent simulation using 3D FE would take hours, if not days, and use many times more memory [START_REF] Pereira | Simulation of acoustic guided wave propagation in cortical bone using a semi-analytical finite element method[END_REF].

The SAFE method can be considered as an accurate tool for simulating the ultrasonic guided wave response. The method has shown results in agreement with the 3D conventional FE simulations, as presented in our previous study [START_REF] Pereira | Simulation of acoustic guided wave propagation in cortical bone using a semi-analytical finite element method[END_REF]. However, SAFE has the advantage of requiring reduced computing resources (e.g. 5.5 h and 2.5 Gb of memory against 4 h and 1.4 Gb of memory for the FE and SAFE methods receptively). Even though the advantages of SAFE in terms of resource requirements are modest, it is worth noticing that once the dispersion curves are calculated, the propagated responses can be computed for different excitation configurations without running the simulations again within seconds at minimal memory cost (<500 Mb).

B. Axial transmission configuration

The axial transmission configuration was modeled using the SAFE procedure described The acoustic excitation was performed using a point source located in a plane denoted z 0 = 0 mm in contact with the upper part of the periosteal region of the bone cross-section.

A 5-cycle Hann windowed toneburst centered at different frequencies between 20-85 kHz was used as the input waveform to perform the excitation. Eight different excitation positions around the circumference of the cortical bone were investigated in this study. The positions, identified with the letters "A" to "H" in Fig. 2, were chosen based on the anatomical accessibility of a probe in a left human forearm. As a limitation of our study, the simulation of the excitability curves was made based on the displacement computed directly on the cortical bone tissue instead of on the overlaying soft tissue (see Fig. 2). The reason for this simplification is due to the higher mesh stability by considering the displacement in the inner elements (cortical tissue) instead of in the outer elements (soft tissue). When considering the displacement on the cortical bone tissue instead of the overlying soft tissue, a negligible change in the excitability curves was observed. However, for modes traveling with most of the energy in the cortical layer (which is the case for the majority of modes investigated in this study), the effect of this assumption on the excitability curves is acceptable. As the wavelengths of the modes of interest are long relative to the soft tissue thickness, a very limited number of modes are traveling with significant energy in the soft tissue. For the modes traveling in the cortical layer, the excitability curves have shown similar shape for both cases so as the relative difference and the rank between modes remained the same. As a consequence, no significant change on the relative variation obtained for the phase velocity and cut-off frequencies in the sensitivity study was observed with or without the soft tissue."

In reception, the out-of-plane displacement of the propagating waves were simulated at from one another by a distance of 3.5 mm, while varied from z = 20 mm to z = 90 mm. The signals obtained were processed using the 2D-FFT [START_REF] Alleyne | A 2-dimensional fourier transform method for the quantitative measurement of lamb modes[END_REF]) in order to calculate the frequency vs. phase-velocity diagram. Figure 3a shows typical received time domain signals obtained for an excitation frequency centered at 45 kHz. Figure 3b shows the frequency vs. phase-velocity diagram obtained after performing the 2D-FFT.

For each diagram, the phase velocity of the low-order mode (named V ph ) was measured by taking the peak of energy at the central frequency of the excitation waveform, as shown in the black square in Fig. 3b. Similarly, the cut-off frequency of the higher-order mode (named F cut-of f ) was measured by taking the peak of energy at a reference velocity equal to 4000 m/s (see black circle in Fig. 3b). Those two features (V ph and F cut-of f ) were used to evaluate the performance of each configuration in term of assessing the variations in the intracortical bone properties. It should be noted that the 2D tri-layer model investigated is this study did not account for the impacts of structural or geometrical irregularities along the axial direction since the waveguide was assumed to be constant and infinite. Such an assumption may change the signal of guided waves propagation in bone tissue. However, this limitation is balanced by the fact that the bone properties are typically measured in the central region of the bone, where the assumption of a constant profile along the axial direction holds because of the long wavelength in the frequency range used in this study.

C. Modeling cortical bone

The geometry of the waveguide was extracted from a slice of a human radius geometry.

The radius geometry was provided by Sawbones (Sawbone, 1997), and was built based on µCT images of the human radius. In this work, the middle third region, which is composed only of cortical bone, was used. Moreover, bone-marrow was added inside the cortical shell, and a 3.5 mm layer of soft-tissue was added outer the cortical shell. The tri-layer crosssection was imported into the COMSOL Multiphysics environment for execution of the SAFE analysis. The mesh was built using triangular elements with a maximum size of 0.5 mm, subject to a constraint of at least 15 elements per wavelength. The bone marrow and soft tissue were modeled as homogeneous viscoelastic fluids with the same properties used by Naili et al. [START_REF] Naili | Influence of viscoelastic and viscous absorption on ultrasonic wave propagation in cortical bone: Application to axial transmission[END_REF] (bulk modulus=2.25e9 Pa, bulk viscosity=1.97 Pa.s and density=1000 kg/m 3 ), while the cortical bone was modeled as a viscoelastic transversally isotropic material. Furthermore, cortical bone was divided into three regions along the radial direction: periosteal (a-b), middle (b-c) and endosteal (d-e) regions (see Fig. 2a). Each TABLE I. Summary of the properties used to model cortical bone. The real part of the stiffness coefficients were adapted from the experimental results obtained by [START_REF] Sansalone | Determination of the heterogeneous anisotropic elastic properties of human femoral bone: From nanoscopic to organ scale[END_REF], while the imaginary parts were adapted from the results presented by [START_REF] Naili | Influence of viscoelastic and viscous absorption on ultrasonic wave propagation in cortical bone: Application to axial transmission[END_REF], using equation 7. et al., 2010). A more detailed description of the distribution of the material properties along the radial direction is provided in the following section.

In the viscoelastic model used in this study, the elastic stiffness constants were written as complex numbers:

C c = C c + η c f i (6)
where the subscripts c=11, 13, 33, 44 and 66, and f corresponds to the frequency. The real part of the stiffness constants is assumed to be independent of the frequency, while the imaginary part results from the use of the linear theory of viscoelasticity, assuming frequency-dependent losses in the waveguide. Because relevant data was lacking, the viscosity coefficients η 11 , η 13 , η 33 , η 44 and η 66 shown in Table 1 were defined based on the viscosity coefficients found in the literature. The coefficients reported in [START_REF] Naili | Influence of viscoelastic and viscous absorption on ultrasonic wave propagation in cortical bone: Application to axial transmission[END_REF] were adjusted by the ratio between the elastic coefficients used in this study and those reported in [START_REF] Naili | Influence of viscoelastic and viscous absorption on ultrasonic wave propagation in cortical bone: Application to axial transmission[END_REF], as given by:

η c = C c C * c .η * c ( 7
)
where the subscripts c=11, 13, 33, 44 and 66, η * c and C * c are respectively the viscosity and elastic coefficients reported by Nailiet al. [START_REF] Naili | Influence of viscoelastic and viscous absorption on ultrasonic wave propagation in cortical bone: Application to axial transmission[END_REF] for the cortical bone.

D. Heterogeneous distribution of bone properties

The heterogeneous nature of the biomechanical properties of cortical bone tissue was adapted from the experimental results obtained by Sansalone et al. [START_REF] Sansalone | Determination of the heterogeneous anisotropic elastic properties of human femoral bone: From nanoscopic to organ scale[END_REF] using 3-D synchrotron micro-computed tomography images. In order to simplify the spatial distribution of each material property, the original values (Fig. 4) were adapted by choosing three different constant values for the periosteal (between a-b), middle (between b-c) and endosteal (between d-e) regions. Moreover, a linear variation was assumed in the region between the middle and periosteal bone (between c-d). For the density, the same approach was applied based on porosity distribution reported in the study. The conversion was made considering a scale where 100% BV/TV (bone volume/ total volume) denotes a bone with a 1850 kg/m 3 density and 0% of BV/TV denotes only bone marrow with a 1000 kg/m 3 density. The constant values, as well as the linear variation of each material property, were chosen to minimize the gap between the original [START_REF] Sansalone | Determination of the heterogeneous anisotropic elastic properties of human femoral bone: From nanoscopic to organ scale[END_REF] and simplified values. Figure 4 shows the original and simplified (healthy and osteoporotic) distributions for the coefficient C 11 .

Two different physiopathological conditions, namely "Healthy" and "Osteoporotic", were defined for the endosteal region (between d*-e* and d-e in Fig. 4, respectively). The "Osteoporotic" condition was defined with the simplified distribution taken from the original data [START_REF] Sansalone | Determination of the heterogeneous anisotropic elastic properties of human femoral bone: From nanoscopic to organ scale[END_REF]. This assumption was made because the data were obtained from an osteoporotic 79-year old patient. The "Healthy" condition was then obtained by: i) arbitrarily increasing the values of all stiffness coefficients and density in the endosteal region by 20%, and ii) by reducing the length of the endosteal region by 20%.

As a limitation, the approach used to model the "Healthy" condition constitutes a simple way of modeling a slightly less degraded condition, which does not necessarily represents the condition associated with a healthy patient. However, since the properties are changing only in the endosteal in our sensitivity study, the magnitude of change between the "Healthy"

and "Osteoporotic" condition can be considered very small, which may be similar to the magnitude of degradation associated to the early-stage osteoporosis. Furthermore, the het- erogeneous nature of cortical bone is based on homogenized local material properties, which are not easy to define due to the multi-scale nature of bone and its dependence on the microstructure at smaller scales. However, such discrepancies tend to be minimized by using low-frequency excitation associated with long wavelengths of the excited modes. Table I shows a summary of the "Healthy" and "Osteoporotic" elastic coefficients and densities used to defined the aforementioned scenarios.

Figures 5a to h show the normalized excitability curves (in the 10-100 kHz range) for an excitation performed at eight different positions on the bone surface (identified with letters "A" to "H" in Fig. 2). The normalized excitability predicts the relative amplitude expected to be observed by a receiver when an emitter at position "A" applies a load in the radial direction (see Fig. 2). The curves show the modal excitability of five modes with similar velocities (m=1, m=2, m=3, m=4, and m=5, marked with a solid line ellipse in Fig. 3).

A notable difference in the modal excitability can be observed between each configuration, which is expected for non-symmetric waveguides such as the bone geometry modeled in this study. For practical applications, a position that can mainly excite a single mode is preferable in order to avoid interferences from other modes.

The modes m=1, m=2, m=3, m=4, and m=5 all present flexural-like mode shapes. For instance, Fig. 6a and 6b show the out-of-plane component of the fundamental flexural tube mode F(1,1) and mode shape of m=1 respectively. Despite the similarities between the displacement fields and the mode order, a direct association to the mode shape of tubular waveguides is not possible due to the non-symmetric nature of the waveguide and its arbitrary geometry

The modes m=1, m=2, and m=3 all show an excitability 10 dB higher than the other modes, as shown in the highlighted regions in Fig. 3a, 3c, and 3h, respectively. Furthermore, according to Figure 7, the modes m=1, m=2, and m=3 have a higher percentage of the total power flow concentrated in the endosteal region as compared to the other modes, which indicates that they are likely to be sensitive to changes in this region. Thus, the excitation at position "A", "C" and "H" were chosen as the most suitable positions to excite a dominant mode; they were therefore further investigated in this paper. Since the attenuation of the modes is expected to be small at very low frequencies (<100 kHz), it was not taken into account to choose the most suitable excitation and detection positions

Effect of physiopathological conditions

Table II shows the percentage variation in the phase velocity (V ph ) and cut-off frequency (F cut-of f ) for the excitation performed at position "A", "C" and "H". The variations was calculated based on the difference computed between the two physiopathological conditions ("Healthy" and "Osteoporotic"). The simulations were performed using input waveforms centered at different frequencies between 20-85 kHz. This choice was made in order to cover a frequency range that has the potential to excite separately one of the five modes with similar velocities on the excitability curves in Fig. 5. The sensitivity of phase velocity (V ph ) was strongly affected by the excitation position.

The highest sensitivity was achieved at position "A", for all frequencies. On the other hand, the cut-off frequency (F cut-of f ) was only slightly affected by the position of excitation,

showing similar sensitivities for all positions. However, the cut-off frequency was strongly affected by the frequency of excitation, showing a significantly higher sensitivity at the first cut-off frequency (20 kHz) as compared to the second (40 kHz) and third (70 kHz). Similarly, an increase in the sensitivity of the phase velocity (V ph ) was observed at lower frequencies, reaching a maximum variation at 45kHz, for all positions.

Figure 7 shows the power flow in the endosteal region (normalized to the total power flow) for modes m=1, m=2, m=3, m=4, and m=5, varying from 20 kHz to 100 kHz. The results show that mode m=1 has a higher power flow in the endosteal region as compared to the other modes, achieving approximately 23% of the total power in this region at 45 kHz.

This explains the higher sensitivity obtained for the excitation performed at position "A", in which the mode m=1 excitability is 15 dB higher than the other modes (see Fig. 5a)

Effect of stiffness coefficients and density

Table III shows the percentage variation computed for the V ph and F cut-of f for each viscoelastic coefficient and density separately. The condition varied from "Healthy" to the "Osteoporotic". The excitation was performed at position "A" using central frequencies ranging between 20-85 kHz.

The sensitivity of V ph showed a positive variation associated with coefficients C 11 , C 33 , C 44 , C 66 , while a negative variation was associated with the coefficient C 13 and the density.

Coefficient C 44 showed the major contribution for most of the frequencies investigated (except at 35 kHz), accounting for approximately 50% of the variations in the degradation of all parameters taken together (identified as "Overall" in Table III). The reduction in sensitivity observed at 35 kHz was found to be associated with the sharp decrease in the power flow of Power flow in the endosteal region normalized by the total power flow for modes m=1, m=2, m=3, m=4, and m=5 from 20 kHz to 100 kHz. mode m=1 seen in Fig. 7 at frequencies below 45kHz. In addition, it may be related to the reduction in excitability seen in mode m=1 at lower frequencies (see the excitability curves in Fig. 5a), which may increase the interference from other modes, and consequently, reduce the sensitivity. For the sensitivity of F cut-of f , a positive variation associated with the coefficient C 66 was observed, accounting for almost 100% of the "Overall" variation. Such a dependence could be used to assess the status of coefficient C 66 separately.

IV. CONCLUSION

In this paper, the SAFE method was used to simulate the propagation of guided waves at low frequency (20-85 kHz) in a bone waveguide in the context of an axial transmission configuration. The SAFE method was shown to be an efficient tool for investigating the effect of intracortical properties on the propagation of guided waves since a number of configurations were simulated using reduced computing resources.

Considering the actual bone geometry, sixteen guided wave modes were found in the 10-100 kHz frequency range (Figure 3b). However, only some of the modes have shown enough excitability and consequently clinical interest. In addition, the excitability and sensitivity of each mode were shown to vary according to the frequency and position of excitation on the bone surface. The phase velocity showed a maximum sensitivity at 45 kHz, which would appear to be associated with the peak of energy (power flow) concentrated in the endosteal region at this frequency. The cut-off frequency was only slightly affected by the position of excitation, but strongly by the frequency of excitation. A maximum sensitivity was achieved for the first cut-off frequency at approximately 20 kHz.

The sensitivity of phase velocity was associated with the variation in the physiopathological conditions of all the coefficients and density, but mainly with C 44 . The sensitivity of the cut-off frequency for its part was mainly associated with variation in the physiopathological conditions of the stiffness coefficient C 66 , allowing the axial shear modulus to be assessed separately.

Thus, by selecting an adequate frequency and position of excitation, the low-frequency axial transmission was shown to be a promising method for assessing intracortical bone properties. The features identified in this study could be used as a metric to compare the similarity between experimental and numerical data. A cost function based on these features could then be implemented into an inversion scheme to retrieve reliable bone properties from experimental data at low frequency. A parameterized bone-like geometry model, instead of a plate or cylinder, could potentially reduce the errors commonly encountered with in-vivo and ex-vivo experiments. For instance, the model could potentially compensate geometrical variations between different patients that are unrelated to the bone quality, such as the outer diameter and the external shape of the bone. The method could, therefore, be applied to identify small changes associated with early-stage osteoporosis or gradual evolution of the bone condition over time. For these extreme cases, since the sensitivity associated to slightly different degraded conditions was found modest, a good basis for comparison (e.g., a baseline method or a large reference database) would have to be implemented in order to guarantee the robustness of the method.

  FIG. 1. Schematic of the equivalent 3D axial transmission configuration used to perform the SAFE

  FIG. 2. Cross-section of the middle of the radius bone used for the SAFE analysis. Different regions were

  FIG. 3. a) Simulated time domain signals obtained at each receiving position. b) Typical phase velocity

  *-e*) Osteoporotic (d-e)*

FIG. 4 .

 4 FIG. 4. Simplified distribution of the elastic coefficients and density varying with the distance from the

  FIG. 5. Normalized modal excitability obtained for the excitation performed at eight different position

FIG. 6 .

 6 FIG. 6.

Frequency

  FIG. 7.

  region was defined with different values for the transversally isotropic stiffness coefficients and density (C 11 , C 13 , C 33 , C 44 , C 66 and ρ), as summarized in Table I. The value of each region was chosen based on a simplification of experimental results reported by Sansalone et al. (Sansalone

					*
	C 11 (Pa)	13.17e9 + 94.95f i 11.46e9 + 82.62f i 4.78e9 + 34.46f i 3.11e9 + 22.42f i
	C 13 (Pa)	6.60e9 + 91.74f i 6.03e9 + 83.81f i 3.32e9 + 46.14f i 2.65e9 + 36.83f i
	C 33 (Pa)	16.04e9 + 109.2f i 14.31e9 + 97.45f i 6.19e9 + 42.15f i 4.17e9 + 28.39f i
	C 44 (Pa)	3.75e9 + 14.70f i 3.15e9 + 12.34f i 0.87e9 + 3.41f i	0.30e9 + 1.17f i
	C 66 (Pa)	3.62e9 + 14.19f i 2.99e9 + 11.72f i 0.78e9 + 3.05f i	0.23e9 + 0.90f i
	Density (kg/m 3 )	1850.00	1850.00	1299.04	1161.30

*The "Healthy" conditions were defined based on an arbitrary gain of 20% on the value of all stiffness coefficients and density. **Adapted from the experimental results obtained by Sansalone et al.

[START_REF] Sansalone | Determination of the heterogeneous anisotropic elastic properties of human femoral bone: From nanoscopic to organ scale[END_REF]

TABLE II .

 II Total variation of the phase velocity V ph and cut-off frequency F cut-of f for the excitations performed at position "A", "C" and "H".

	Phase velocity (V ph ) variation (%)	Position A Position C Position H
	35 kHz	5.59	4.91	4.21
	45 kHz	7.98	6.01	4.23
	55 kHz	7.36	6.23	4.03
	65 kHz	6.72	5.44	2.18
	75 kHz	6.21	4.99	2.50
	85 kHz	5.99	4.48	2.07
	Cut-off frequency (F cut-of f ) variation (%) Position A Position C Position H
	First cut-off frequency (20 kHz)	17.37	16.31	16.31
	Second cut-off frequency (40 kHz)	9.42	12.15	11.23
	Third cut-off frequency (70 kHz)	5.45	6.99	5.38
	*The variations were calculated based on the difference computed between the "Healthy" and the
	"Osteoporotic" condition.

TABLE III .

 III Percentage variation computed for V ph and F cut-of f for each viscoelastic coefficient and density separately.

					Position A
	Phase velocity (V ph ) variation (%)	C 11 C 13 C 33 C 44 C 66 Density Overall
	35 kHz	0.32 -1.42 0.8 1.30 1.63 -0.79	5.59
	45 kHz	0.37 -1.45 1.24 3.30 2.01 -0.611 7.98
	55 kHz	0.31 -2.23 1.47 3.21 0.31 -0.72	7.36
	65 kHz	0.36 -2.02 1.73 3.25 0.36 -0.62	6.72
	75 kHz	0.16 -1.73 1.54 3.25 0.32 -0.71	6.21
	85 kHz	0.08 -1.14 1.20 3.21 0.42 -0.67	5.99
	Cut-off frequency (F cut-of f ) variation (%)
	First cut-off frequency (20 kHz)	0	0	0	0 22.06 -2.34 17.37
	Second cut-off frequency (40 kHz)	0	0 0.90 0.90 9.42 -0.90	9.42
	Third cut-off frequency (70 kHz)	0.94 0.52 0.41 0.94 3.98 -1.05	5.45

ACKNOWLEDGMENTS

This work was developed in collaboration with the Centre de recherche l'Hôpital du Sacré-Coeur de Montréal and École de technologie supérieure. The authors acknowledge the grant from the Fonds d'Internationalisation de la Recherche (FIR) of École de technologie supérieure.


Philippe Monmousseau
email: philippe.monmousseau@enac.fr

Daniel Delahaye

Aude Marzuoli
email: amarzuoli3@gatech.edu

Eric Feron
email: feron@gatech.edu

Door-to-door Travel Time Analysis from Paris to London and Amsterdam using Uber Data

Keywords: Door-to-door travel times, Big Data, Multi-modal travel

SESAR Flightpath 2050 sets some ambitious goals for air travel, specifically regarding door-to-door travel times. Using Uber data, a reliable estimation of door-to-door travel times is possible which then enables a comparison of different travel modes such as rail and air. This model of the full door-to-door travel time for multi modal trips can also be used to evaluate where Europe stands with respect to SESAR's goals and how they should try to address the related difficulties.

I. INTRODUCTION

Seamless door-to-door travel and data sharing was deemed as needed by the European Commission's 2011 White Paper [START_REF] Commission | White Paper: Roadmap to a Single European Transport Area -Towards a Competitive and Resource Efficient Transport System[END_REF]. Data sharing was already a main focus in the early 2000s and led at an air system level to the creation of the architecture SWIM -System Wide Information Management [START_REF] Meserole | What is System Wide Information Management (SWIM)?[END_REF]. SESAR Flightpath 2050 [START_REF] Commission | Flightpath 2050: Europe's Vision for Aviation ; Maintaining Global Leadership and Serving Society's Needs[END_REF] aims to take a more passenger-centric approach and sets some ambitious goals, including some that are not measurable yet due to lack of available data. Regarding door-to-door travel times, it aims at having 90% of travellers within Europe being able to complete their door-to-door journey within 4 hours. The shift from flight-centric information to passenger-centric metrics was first explored by Cook et al. [START_REF] Cook | Passenger-Oriented Enhanced Metrics[END_REF] within the project POEM -Passenger Oriented Enhanced Metrics.

In the United States, Marzuoli et al. [START_REF] Marzuoli | Implementing and validating air passenger-centric metrics using mobile phone data[END_REF] presented a method to detect domestic air passengers on a nationwide scale using mobile phone data, enabling a per leg analysis of the full door-to-door trip though the main focus was on passengers' behavior at airports. The passengers' experience in airports under major perturbations using this method and additional data from social media was further studied in [START_REF] Marzuoli | Passenger-centric metrics for Air Transportation leveraging mobile phone and Twitter data[END_REF]. In Europe, within the BigData4ATM project 1 , García-Albertos et al. [START_REF] García-Albertos | Understanding Doorto-Door Travel Times from Opportunistically Collected Mobile Phone Records[END_REF] presented a methodology for measuring the door-to-door travel time using mobile phone data and applied it between two Spanish cities, Madrid and Barcelona. Mobile phone data is however proprietary data and difficult to access for research.

Sun et al. [START_REF] Sun | Competitiveness of on-demand air taxis regarding door-to-door travel time: A race through Europe[END_REF] implemented a door-to-door minimum travel time estimation based on open source maps and datasets in 1 www.bigdata4atm.eu order to study the possible competitiveness of air taxis. This study is also based on readily available online data but aims in creating a method to measure the actual average door-todoor travel time once the trips are over enabling an analysis and comparison of the different modes. This method is here applied to two different European trips: from Paris to London and from Paris to Amsterdam.

This paper is organized as follows: Section II presents the model and data used to evaluate the full door-to-door journey time. Section III showcases some use cases of this model and Section IV concludes this paper and discusses further research directions.

II. MODEL

Similarly to [START_REF] García-Albertos | Understanding Doorto-Door Travel Times from Opportunistically Collected Mobile Phone Records[END_REF] and [START_REF] Sun | Competitiveness of on-demand air taxis regarding door-to-door travel time: A race through Europe[END_REF], the full door-to-door travel time can be decomposed into the five following times:

T = t to + t dep + t in + t arr + t from (1)
where • t to is the time needed to go from the start of the journey to the selected modal station (e.g. airport or train station) • t dep is the time spent waiting and going through security processes at the departure station • t in is the time actually spent in the main mode's vehicle • t arr is the time spent going through security processes at the arrival station • t from is the time needed to go from the arrival modal station to the end of the journey The measurement and estimation of these different times are described in the upcoming subsections. This study limits its scope to the two following modes: by air and by rail. The rail mode is then further separated between the Eurostar and the other rail services.

A. Time to and from the modal stations

Uber, a ride-sharing service implanted in major urban areas on six continents, recently released anonymized travel times data for certain of these urban areas. This data consist of the average travel time between zones (e.g. census tracts) within the serviced area from all Uber rides aggregated over each considered day. Depending on the availability of data, five additional different periods are considered:

• Early Morning: from midnight to 7am • Morning: from 7am to 10am • Midday: from 10am to 4pm • Afternoon: from 4pm to 7pm • Evening: from 7pm to midnight Before this data release, several studies were already conducted on the impact of Uber in urban transit. Hall et al. [START_REF] Hall | Is Uber a substitute or complement for public transit?[END_REF] studied the impact of Uber on the use of public transit system based on Uber's entry date in different cities and focused on whether Uber complemented or substituted public transit. Similarly Li et al. [START_REF] Li | Do On-demand Ride-sharing Services Affect Traffic Congestion? Evidence from Uber Entry[END_REF] concluded that on average Uber tends to decrease congestion in US urban areas where it was introduced. Wang and Mu [START_REF] Wang | Spatial disparities of Uber accessibility: An exploratory analysis in Atlanta, USA[END_REF] studied Uber's accessibility in Atlanta, US by using the average wait time for a ride as a proxy and concluded that the use of Uber was not associated to a specific social category. Uber rides being part of the road traffic flow, this study considers Uber's travel times as accurate proxies of the actual travel time by road (whether by car or by bus). The scope of this paper limits itself to road access and egress to the modal stations considered, though subway alternatives should be considered by using schedules and real time data for a more exhaustive analysis of these legs of multi-modal trips.

Data was gathered from Uber's Movement API2 over the period of January 1st 2018 to March 31st 2018. Paris was divided into the IRIS zones used by INSEE for census, London into the Lower Super Outer Area (LSOA) division defined by the Greater London Authority, and Amsterdam into its official districts called wijk.

B. Time at the modal stations

For this study, the time spent at each modal station is considered to depend only on the mode and whether it is at boarding or unboarding process. Those times are based on the train stations websites for rail stations, and on the study of [START_REF] Marzuoli | Implementing and validating air passenger-centric metrics using mobile phone data[END_REF] for the time spent at airports. They are summarized in the following table: The difference between time spent at departure at a normal rail station and at the Eurostar station comes from the necessity of passport checking since the United Kingdom is not part of the Schengen passport-free area. Better estimates could be obtained using data gathered from GPS or mobile phone sources as well as WiFi beacons within airports.

C. Time in main mode

For this study, it was assumed that the modes were on time and based on a weekly schedule. Future studies should consider using actual flight or train times, taking into account delays and perturbations. The Eurostar schedule was taken from its website 3 and is shown in Table II. The rail schedule between Paris and Amsterdam was based on a week's schedule from the month of July 2019 gathered using an online train booking service 4 and is summarized in Table III. Similarly, the flight schedules from Paris Charles de Gaulle airport (CDG) to London Heathrow airport (LHR) and to Amsterdam Airport Schiphol (AMS) were extrapolated from the flight schedule of the third week of January 2019 gathered from CDG's website 5 . Both schedules can be found in Tables IV & V. Using the updated schedules with delays and cancellations would be a promising next step and easily integrated to this model.

D. Full door-to-door time

For this study, it was assumed that the main mode of transport would leave on time and that travellers planned their departure time to arrive at the departure station exactly t dep minutes before the scheduled departure time. This assumption is necessary to determine which period of the day to consider when retrieving the Uber average time between the initial zone and the departure station to determine t to . For days and zones where only daily aggregates were available in the Uber data, the daily aggregated times were used for each period of the day as a proxy. The same process was implemented to select the value of t from .

III. APPLICATIONS

A. Modal comparison

Using this model and the available data makes it possible to analyze at what time of the day it is usually best to start the journey from Paris. This analysis enables an easy hourly modal comparison as well. As an example, this study considers the case of a tourist visiting the Eiffel Tower in Paris and having to go back to London or Amsterdam directly from there. Travel times for each scheduled flight and each scheduled train were extracted for the period of January to March 2018.

1) From Paris to London: Fig. 1 shows the hourly boxplots of the full door-to-door travel time for a Londoner using either the Eurostar train or flying from CDG to LHR. The four hour goal is marked in red. A first observation is that this time goal is for the moment reached only by taking the train, even though most of London's LSOAs are not reachable within that time limit. Regarding the train, there are only two hours in the day that fare poorly compared to the otherwise constant behavior of this mode: the 9am rush and at 4pm, just before the closing of businesses. Fig. 2 shows the boxplots per day of the week of these full door-to-door travel time along with the four hour goal mark. Regarding this time limit, it is achieved every day of the week via train for some areas of London. These boxplots also indicate that the best days to travel to London from the center of Paris are not the same depending on the chosen mode. When taking the plane, it is best to fly on Saturdays and Sundays, while it is best to take the train on Mondays, Tuesdays and Saturdays.

Using this model, it also possible to determine which zones in the destination city can be most easily reached depending on the chosen mode. Fig. 3 shows which zones were reached using Uber rides from the train station on January 8, 2018 at two darker greens are the zones reachable in four hours or less in the Eurostar case (Fig. 3), which corresponds to the center of London. By plane, the fastest time is 4h30 for the surrounding of the airport, while reaching the center of London takes around five hours.

2) From Paris to Amsterdam: Fig. 5 shows the hourly boxplots of the full door-to-door travel time for a Amsterdamer using either the train or flying from CDG to AMS. A first observation is that the four hour limit is not reached by either mode. Taking the train brings you closer to this goal than taking the plane in the morning, except at rush hour 9am. In the afternoon it is fastest to take the plane, except for those who prefer night trains in order to arrive early in the morning of the next day. the train shown in these boxplots: While the full door-to-door travel time distribution is nearly constant by plane, there is a major difference between traveling by rail during weekdays and on the weekend: the median traveling time on weekdays is higher by around ninety minutes than on weekends. Fig. 7 shows which zones were reached using Uber rides from the train station on January 11, 2018 at 17:50, these zones being colored based on the full door-to-door travel time needed from the Eiffel Tower. Fig. 8 shows the same map for air travelers on January 11, 2018 at 08:00 but with a different color scaling. For these dates, taking the plane is faster than the train for every considered wijk.

B. SESAR Flightpath 2050 objectives

With the Flightpath 2050 in mind, this model can also be used to assess on which phase of air travels improvements should be made, assuming the 4 hour door-to-door travel time limit is to be achieved by means of air transportation. The following analysis could be similarly conducted to assess rail transportation system. Fig. 9 shows the boxplots of the time spent in the different segments using air and rail in order to reach to London assuming the traveler is leaving from the Eiffel Tower area. It is to be noted that the time spent leaving the airport has a greater variation than the time going to the airport due to the fact that the starting zone is fixed (the Eiffel Tower area) whereas all zones reached by Uber rides from the arriving airport are considered. By rail, the traveler spends the most time in the train, whereas an air traveler spends more time waiting at the departure airport than in flight and can spend more time reaching its destination in London than in any other segment. Improving airports processes at departure and thus decreasing the time spent within the departure airport could be a first step in reaching the 4 hour goal. Fig. 10 shows the boxplots of the time spent in the different segments using air and rail in order to reach to Amsterdam assuming the traveler is leaving from the Eiffel Tower area. Unlike the trip to London, the air traveler spends less time leaving the airport than for any other phase. The 90 minute wait time at the departure airport is here again the most consuming. Regarding rail travels, the time spent in the train is dominant compared to any other travel phase and needs to be addressed if the 4 hour limit is to be applied to rail transit as well. This model has assumed fixed wait times at airports based on a recent study [START_REF] Marzuoli | Implementing and validating air passenger-centric metrics using mobile phone data[END_REF], yet these times can be modified in order to assess the impact of improving the airport's processing of passengers. The following examples assumes that the wait time at arrival drops at an average of thirty minutes and investigates the drop necessary in the processing time at the departure airport. Fig. 11 shows the hourly boxplots of the full door-todoor travel time for a Londoner flying from CDG to LHR using three different wait times at the departure airport: 30, 60 and 90 minutes. Clearly dropping only the processing time at arrival is not enough for this trip. Reducing the processing time at the departure airport to 30 or 60 minutes enables to reach some zones in London in less than 4 hours, but never the totality. And for five rush hours (8am, 10am, noon, 5pm and 7pm), these drops are not enough to reach any of London's LSOA in less than 4 hours. 12 shows the hourly boxplots of the full door-to-door travel time for a Amsterdamer flying from CDG to AMS using the same three different wait times at the departure airport. As for London, only reducing the arrival processing time is not enough to reach the 4 hour time limit. However, reducing the processing time at the departure airport to an hour is enough to reach most of Amsterdam's wijken with the exception of four rush hours (1pm, 3pm, 5pm and 7pm). Decreasing further the processing time to half an hour ensures that the 4 hour limit is reached for the quasi majority of wijken, with the exception of the same four rush hours. This first study indicates that reducing the processing time at airports is a necessary step to achieve the Flightpath 2050 goal, however it is not sufficient: access to and from the airports has to be improved since rush hours can jeopardise the overall achievement of this 4 hour limit. It is to be noted that the effect of these rush hours could be mitigated by the use of other types of access modes such as the subway. This paper leveraged Uber's recently released data in order to create a model of the full door-to-door travel time between two pairs of European cities using either the train or the plane. It could however be implemented for any world city pairs with available ride-sharing or taxi data. This model enables the comparison of the different traveling modes between two cities on at least two different levels: a first approach is to compare the time spent in the different phases of the trip, while a second approach assumes the existence of a time limit and compares the reach of each considered mode. This model also has two additional benefits with respect to SESAR's Flightpath 2050: it can be used both to evaluate where Europe stands regarding its objectives as well as which directions to consider in order to reach these goals.

Further studies should consider using actual schedules of the different modes when evaluating the full door-to-door travel time, as well as alternative modes to reach the airport or train stations such as the subway. Additionally, knowing the daily proportion of travellers using the different approaches would enable a single daily evaluation of the full door-to-door travel time. A possible method to determine this proportion would be by using aggregated information from GPS or mobile phone sources.

Figure 1 :

 1 Figure 1: Hourly boxplots of door-to-door travel times to London starting from the Eiffel Tower in Paris. The top figure shows the boxplot related to air travel, the bottom figure the ones related to train travel.

Figure 2 :

 2 Figure 2: Day-of-week boxplots of door-to-door travel times to London starting from the Eiffel Tower in Paris. The top figure shows the boxplot related to air travel, the bottom figure the ones related to train travel.

Figure 3 :

 3 Figure 3: Door-to-door travel times to London starting from the Eiffel Tower in Paris using the Eurostar train

Figure 4 :

 4 Figure 4: Door-to-door travel times to London starting from the Eiffel Tower in Paris by plane via CDG and LHR

Figure 5 :

 5 Figure 5: Hourly boxplots of door-to-door travel times to Amsterdam starting from the Eiffel Tower in Paris. The top figure shows the boxplot related to air travel, the bottom figure the ones related to train travel.

Fig. 6

 6 Fig.6shows the boxplots per day of the week of these full door-to-door travel time along with the four hour goal mark. There is a striking contrast between taking the plane and taking

Figure 6 :

 6 Figure 6: Day-of-week boxplots of door-to-door travel times to Amsterdam starting from the Eiffel Tower in Paris. The top figure shows the boxplot related to air travel, the bottom figure the ones related to train travel.

Figure 7 :

 7 Figure 7: Door-to-door travel times to Amsterdam starting from the Eiffel Tower in Paris using the train

Figure 8 :

 8 Figure 8: Door-to-door travel times to Amsterdam starting from the Eiffel Tower in Paris by plane via CDG and AMS

Figure 9 :

 9 Figure 9: Travel time boxplots for the different phases of travel to London starting from the Eiffel Tower in Paris. The top figure shows the boxplot related to air travel, the bottom figure the ones related to train travel.

Figure 10 :

 10 Figure 10: Travel time boxplots for the different phases of travel to Amsterdam starting from the Eiffel Tower in Paris. The top figure shows the boxplot related to air travel, the bottom figure the ones related to train travel.

Figure 11 :

 11 Figure 11: Hourly boxplots of door-to-door air travel times to London starting from the Eiffel Tower in Paris for varying wait times at the departing airport

Figure 12 :

 12 Figure 12: Hourly boxplots of door-to-door air travel times to Amsterdam starting from the Eiffel Tower in Paris for varying wait times at the departing airport

TABLE I :

 I TIME SPENT AT MODAL STATIONS

	Mode	Time at departure	Time at arrival
	Air	90 min	45 min
	Rail	15 min	10 min
	Eurostar	45 min	10 min

TABLE II :

 II EUROSTAR WEEKLY SCHEDULE (LOCAL TIME) FROM PARIS TO LONDON

	Mon	Tue	Wed	Thu	Fri	Sat	Sun	Paris	Lon.
	x	x	x	x	x	x		07:13	08:32
	x	x	x	x	x			07:43	09:00
						x	x	08:13	09:30
	x					x		08:43	10:00
	x	x	x	x	x	x	x	09:04	10:39
		x	x	x	x	x	x	10:13	11:30
	x							10:13	11:43
	x				x		x	10:43	12:09
	x							11:13	12:39
		x	x	x	x	x		11:13	12:47
							x	11:28	13:09
	x	x	x	x	x	x	x	12:13	13:30
	x	x	x	x	x	x	x	13:13	14:39
							x	13:28	15:00
						x		14:13	15:39
	x			x	x		x	14:43	16:02
	x	x	x	x	x	x	x	15:13	16:39
							x	16:43	18:02
	x	x	x	x	x	x		17:13	18:32
							x	17:43	19:00
	x	x	x	x	x		x	18:13	19:39
		x	x	x				18:43	20:02
	x	x	x	x	x	x	x	19:13	20:39
					x		x	19:43	21:09
	x	x	x	x		x		20:13	21:49
					x		x	20:43	22:00
	x	x	x	x	x		x	21:13	22:39

TABLE III :

 III THALYS WEEKLY SCHEDULE (LOCAL TIME) FROM PARIS TO AMSTERDAM

	Mon	Tue	Wed	Thu	Fri	Sat	Sun	Paris	Ams.
	x	x	x	x	x			06:13	09:44
	x	x	x	x	x	x		07:25	10:44
	x	x	x	x	x	x	x	08:25	11:44
	x	x	x	x	x	x	x	10:25	13:44
	x	x	x	x	x			11:25	14:44
	x	x	x	x	x			11:55	15:44
	x	x	x	x	x			12:46	17:35
	x	x	x	x	x			12:52	18:35
	x	x	x	x	x			14:25	19:35
	x	x	x	x	x	x		15:55	20:35
	x	x	x	x	x	x	x	17:55	22:35
	x	x	x	x	x	x		18:25	23:38
	x	x	x	x	x			20:55	08:35
						x	x	12:25	15:44
						x	x	12:55	17:35
						x	x	14:25	17:44
						x		19:25	23:41
							x	16:25	19:44

TABLE IV :

 IV SIMULATED WEEKLY SCHEDULE FOR FLIGHTS BETWEEN CDG AND AMS (LOCAL TIME)

	Mon	Tue	Wed	Thu	Fri	Sat	Sun	CDG	AMS
	x	x	x	x	x	x	x	07:10	08:30
	x	x	x	x	x		x	08:05	09:25
						x		08:15	09:35
		x		x				08:55	10:15
	x		x		x	x	x	09:05	10:25
	x	x	x	x	x			09:35	10:55
	x	x	x	x	x	x	x	10:20	11:40
	x	x	x	x	x	x	x	11:45	13:05
	x	x	x	x	x	x	x	12:40	14:00
	x	x	x	x	x	x	x	14:30	15:50
	x	x	x	x	x		x	16:55	18:15
	x	x	x	x	x	x	x	18:10	19:30
	x	x	x	x	x	x	x	18:40	20:00
	x	x	x	x	x	x	x	20:25	21:45

TABLE V :

 V SIMULATED WEEKLY SCHEDULE FOR FLIGHTS BETWEEN CDG AND LHR (LOCAL TIME)

	Mon	Tue	Wed	Thu	Fri	Sat	Sun	CDG	LHR
					x	x	x	07:15	07:40
	x	x	x	x	x	x	x	07:20	07:45
		x						07:30	07:55
	x		x	x				07:35	08:00
	x	x	x	x	x	x	x	10:00	10:25
				x	x			10:20	10:45
	x					x	x	10:25	10:50
		x	x					10:35	11:00
	x		x		x	x	x	11:45	12:10
		x		x				11:50	12:15
	x	x	x	x	x	x	x	13:20	13:45
		x		x			x	14:55	15:20
	x		x					15:05	15:30
					x			15:20	15:45
	x	x	x	x	x	x		16:10	16:35
	x	x	x	x	x	x	x	16:15	16:40
	x	x	x	x	x	x	x	18:05	18:30
	x	x	x	x	x		x	20:35	21:00
	x	x	x	x	x	x	x	21:20	21:45

movement.uber.com

www.eurostar.com

www.thetrainline.com

https://www.parisaeroport.fr

ACKNOWLEDGMENT

The authors would like to thank Nikunj Oza from NASA-Ames, the BDAI team from Verizon Media in Sunnyvale, California, as well as the French government for their financial support.