Ankit Jain
email: ankit.jain@telecom-bretagne.eu

Pascal Pagani

Rolland Fleury

Michel M Ney

Patrice Pajusco

Efficient time domain HF

Keywords: HF propagation, geolocation, ionosphere, multiple recivers, TDoA. I

teaching and research institutions in France or abroad, or from public or private research centers.

INTRODUCTION

HF radio waves can travel long distances using skywaves and therefore have their importance in military and civil applications. In such contexts, it is interesting to locate the position of the transmitter (Tx) by using the different types of information acquired from the signals received at the receiving station. Different methods are proposed to find the Direction of Arrival (DoA) of the signals by deploying a large array of co-located antennas [START_REF] Gething | HF direction finding[END_REF]. Using time domain HF geolocation techniques such as Time of Arrival (ToA) and TDoA, signals received at different non co-located single antenna receivers (Rx's) are sufficient to locate the position of the Tx [START_REF] Pagani | Time domain HF geolocation: Experimental measurements and preliminary results[END_REF]. The ToA method requires a minimum of three Rx's to locate the position of the Tx [START_REF] Bensky | Wireless positioning technologies and applications[END_REF]. Using the propagation path lengths obtained from the three Rx's, the Tx location can be retrieved using trilateration. If the transmission time of the signal is unknown, the Tx position can be retrieved by multilateration using TDoA at different Rx's. The TDoA method requires a minimum of four Rx's to locate the position of the Tx [START_REF] Torrieri | Statistical theory of passive location systems[END_REF].

A long range HF radio transmission uses skywaves where the signal is reflected from the ionosphere. The signal reflected once from the ionosphere back to the earth can be referred as single hop mode. Using a single hop mode, transmissions using skywaves are possible up to distances of 3000 km at low elevation angles. In order to obtain the HF Tx location, it is necessary to obtain the group path lengths between the Tx and the Rx's which largely depend upon the electron density profile of the ionosphere. The skywaves are reflected back to the earth due to the presence of electrons in the ionosphere. The ionization in the ionosphere varies a lot with the solar cycle, geographical location, seasonally and diurnally. These ionospheric variations are the main issues in HF Tx geolocation through the skywaves. HF Tx geolocation can be further improved by the use of prediction models since it reduces time delay errors introduced by the ionosphere [START_REF] Rao | Prediction of differential time-delay errors in HF hyperbolic position-fixing systems[END_REF].

In our previous study of ToA geolocation [START_REF] Pagani | A study of HF transmitter geolocation through single-hop ionospheric propagation[END_REF], the Tx location was obtained without the knowledge of reflection height and the minimum number of 3 required Rx's. The geolocation error was in the range of about 50-100 km depending on whether the transmitter was very close or far away from the Rx's. Another similar variant of HF geolocation was explained later in [START_REF] Wang | Passive TDOA and DOA based HF geolocation without ionosphere information[END_REF] based on TDoA and DoA methods which use multipath and achieve good geolocation performance without any ionospheric information. However, there is a requirement of antenna array on each reception site; accurate TDoA measurements are difficult to obtain in realistic scenario.

This paper studies, the effect of increasing the number of Rx's on the performance of HF Tx geolocation in the radius of 2000 km. HF communication in the span of 2000 km is possible using a single hop mode. Further, to take into account the variable density of electrons in the ionosphere, a QP model of the ionosphere [START_REF] Croft | Exact ray calculations in a quasiparabolic ionosphere with no magnetic field[END_REF] is used. This model provides realistic group path lengths from which corresponding time delays are evaluated. Section II in the paper explains different aspects of ionospheric propagation. Section III explains the geolocation algorithm for TDoA along with its mathematical formulations. Simulation results are presented in Section IV in order to assess the Tx geolocation performance with an increasing number of receivers. Lastly conclusions are drawn in Section V.

II. IONOSPHERIC PROPAGATION:PHYSICS AND MODELING

A. Structure of the Ionosphere

The ionosphere [START_REF] Davies | Ionospheric radio[END_REF] is a region extending from a height of about 50 km to 2000 km over the surface of the Earth. The neutral atoms present in this region are ionized by extreme ultraviolet (EUV) light from the sun. It results in the formation of free electrons and positively charged ions. These free electrons formed by photo-ionization are mainly responsible for the reflection of HF signal from different layers of the ionosphere facilitating long distance transmission.

Due to the presence of different gases at different altitudes in the ionosphere, there exist ionospheric layers. They are known as the D, E, F 1 and F 2 layer in the order of increasing altitude. The F 2 layer is sometimes just referred as the F layer. Due to the opposing behavior between the density of the neutral atmosphere and the intensity of EUV radiation from the sun, the electron density of a given ionospheric layer is a convex function of the altitude, first increasing from zero to a maximum electron density, then again decreasing down to zero. The electron density is the least in the D layer and the maximum in the F layer.

The D layer exists at an altitude of about 50-100 km and is mainly responsible for the attenuation of the radio waves. The absorption or the attenuation of the signal reduces as the frequency of the signal is increased. The E layer exists from 100 to 150 km and can reflect HF signals. In the case of HF radio propagation, the F layer is the most important region. The height range for the F region is the highest which allows large communication paths. Furthermore, the electron density being the highest in the F region, it allows the reflection of the highest frequencies in the HF range.

B. Analysis of the QP model

To take account for the electron density profile for the radio signals propagating through the ionosphere, the QP model can be used. Developed by Croft and Hoogasian [START_REF] Croft | Exact ray calculations in a quasiparabolic ionosphere with no magnetic field[END_REF], it is obtained by a slight modification of the parabolic model used to define the electron density profile of the ionosphere. The QP model provides analytic expressions describing ray path parameters of oblique propagation considering the ionospheric medium. Using the QP model, the ground distance and group path distance can be evaluated between different Tx-Rx's pair for a specific elevation angle (β) from the Tx and operating frequency (𝑓𝑓) by the formulae provided in [START_REF] Croft | Exact ray calculations in a quasiparabolic ionosphere with no magnetic field[END_REF]. The group path is obtained by multiplying the signal transmission time with the speed of light.

The electron density profile of a layer in the QP model is defined as follows:

else , 0 , 1 2 2              - ≤ ≤                       - - = m b b m b b m m m e y r r r r r r r y r r N N (1)
where 𝑁𝑁 𝑒𝑒 is the electron density having a maximum value 𝑁𝑁 𝑚𝑚 at a maximum distance of 𝑟𝑟 𝑚𝑚 , 𝑟𝑟 is the radial distance from the Earth's centre being equal to the sum of the Earth's radius and the height at which reflection occurs in the ionosphere, 𝑟𝑟 𝑏𝑏 is the value of 𝑟𝑟 at the base of the layer and 𝑦𝑦 𝑚𝑚 is the layer's semi-thickness [START_REF] Croft | Exact ray calculations in a quasiparabolic ionosphere with no magnetic field[END_REF].

The QP model of the ionosphere is a monolayer model providing ray path parameters considering reflection only from a specific layer in the ionosphere. It is a limitation of the QP model. Fig. 1 represents the ground distance (D) plotted using QP model at different values of 𝐹𝐹, which is the ratio of the operating frequency (𝑓𝑓) of the source of HF radio waves to the critical frequency (𝑓𝑓 𝑐𝑐 ) of the QP layer [START_REF] Croft | Exact ray calculations in a quasiparabolic ionosphere with no magnetic field[END_REF]. The maximum ground distance at different values of 𝐹𝐹 can be calculated at the elevation angle β = 0 as seen in Fig. 1. Thus, 𝐷𝐷(β = 0) is the maximum distance at which a Tx-Rx link could be established using single hop propagation through an ionospheric layer following the QP model.

The skip distance arises when 𝑓𝑓 𝑐𝑐 of an ionospheric layer is greater than 𝑓𝑓 of the radio signal [START_REF] Sales | High Frequency (HF) Radiowave Propagation[END_REF]. Hence, skip distance occurs only when the value of 𝐹𝐹 > 1. Skip distance can be found by calculating the minimum value of the ground distance curve for any specific value of 𝐹𝐹 > 1. Furthermore, it can also be found by computing the ground distance at which the slope of the curve is infinite for a plot of ground distance versus elevation angle. It is not possible to communicate using HF radio signals in the skip range using skywaves. The ground distance starts increasing again after reaching the skip distance for increasing elevation angle. It is due to high rays, also known as Pedersen rays, whereas the ground distance at which communication can be established before the skip distance range is due to low rays. Skip distance is a point where both the high and the low rays of a HF signal coincide.

III. HF GEOLOCATION ALGORITHM

A single hop transmission can be assumed for geolocation of HF transmitters in the radius of 2000 km. The signal reaches the destination via a single reflection from the ionospheric layer. Fig. 2 represents the HF geolocation geometry considering a uniform, single layer ionosphere. The transmitter is located at point Tx and the receivers are located at point Rx1 and Rx2. As per the Breit and Tuve's theorem [START_REF] Davies | Ionospheric radio[END_REF], the group path between Rx1-Tx sketched by the red curve is equal to the triangular path Rx1-I1-Tx and the same applies for Rx2-Tx.

Due to the varying refractive index in the ionosphere, the signals propagate along a curved path. The apogee height of the group path length is denoted by ℎ 𝑎𝑎 . An equivalent height ℎ 𝑒𝑒𝑒𝑒𝑒𝑒 corresponding to the reflection height of the path length with the same group delay can be defined assuming the signals are reflected via specular reflection on a flat ionosphere. The ionosphere is assumed to be uniform, so the reflection height for all considered paths is similar for signals having similar elevation angles. The Earth is considered to be flat by neglecting its curvature.

Due to electromagnetic image theory, the length of the propagation path between Rx1-Tx is equal to Rx1-I1-Tx'; similarly the propagation paths between Rx2-Tx is equal to Rx2-I2-Tx'.The geolocation process consists of finding the coordinates of Tx' which is situated at a height of 2ℎ 𝑒𝑒𝑒𝑒𝑒𝑒 vertically above the location of the Tx on a flat Earth. The coordinates of Tx' can be found by multilateration in the case of TDoA (intersection of three hyperboloids) as explained in Section I.

The geometrical problem of identifying the intersection of three spheres or three hyperboloids requires solving systems of quadratic equations which cannot be solved analytically. A mathematical approach to obtain the location of Tx' is explained further for the TDoA case using [START_REF] Bakhoum | Closed-form solution of hyperbolic geolocation equations[END_REF] as a reference.

A. Mathematical approach of TDoA

Let the receivers (Rxi) be located at position 𝐏𝐏 𝐢𝐢 = (𝑥𝑥 𝑖𝑖 , 𝑦𝑦 𝑖𝑖 , 𝑧𝑧 𝑖𝑖 ) where 𝑖𝑖 = 1,2,3. . 𝑛𝑛; let the transmitter (Tx) on the flat Earth surface be located at position 𝐏𝐏 𝟎𝟎 = (𝑥𝑥 0 , 𝑦𝑦 0 , 𝑧𝑧 0 ). The coordinates of point Tx' is located at 𝐏𝐏 𝟎𝟎 ′ = (𝑥𝑥 0 , 𝑦𝑦 0 , 𝑧𝑧 0 ′ ) where 𝑧𝑧 0 ′ is given by (2):

𝑧𝑧 0 ′ = 𝑧𝑧 0 + 2ℎ 𝑒𝑒𝑒𝑒𝑒𝑒 (2)
The propagation path length 𝐷𝐷 𝑖𝑖 between a Tx-Rx pair via one hop mode is given by (3):

𝐷𝐷 𝑖𝑖 = 𝑐𝑐(𝑡𝑡 0 -𝑡𝑡 𝑖𝑖 ) = |𝐏𝐏 𝟎𝟎 ′ -𝐏𝐏 𝐢𝐢 | (3)
where 𝑡𝑡 0 is the signal transmit time, 𝑡𝑡 𝑖𝑖 is the signal arrival time at the i th receiver, 𝑐𝑐 is the speed of light, 𝐷𝐷 𝑖𝑖 is the propagation path length between the transmitter and i th receiver and | | represents the norm. Squaring (3), one obtains an equation which is given by:

|𝐏𝐏 𝟎𝟎 ′ -𝐏𝐏 𝐢𝐢 | 2 = 𝑐𝑐 2 (𝑡𝑡 0 -𝑡𝑡 𝑖𝑖 ) 2 (4)
In the case of TDoA, the transmission time 𝑡𝑡 0 is unknown. Computing the difference between (3) for 𝑖𝑖 = 1 (𝑅𝑅𝑥𝑥 1 ) and ( 3) for 𝑖𝑖 = 2 (𝑅𝑅𝑥𝑥 2 ) ; after further simplification one obtains (5) which is independent of 𝑡𝑡 0 and given by:

-2(𝐏𝐏 𝟐𝟐 -𝐏𝐏 𝟏𝟏 ) ⊤. 𝐏𝐏 𝟎𝟎 ′ 𝑡𝑡 2 -𝑡𝑡 1 + 2𝑐𝑐 2 𝑡𝑡 0 = (𝑡𝑡 1 + 𝑡𝑡 2 )𝑐𝑐 2 - �|𝐏𝐏 𝟐𝟐 | 2 -|𝐏𝐏 𝟏𝟏 | 2 � 𝑡𝑡 2 -𝑡𝑡 1 (5)
where ⊤ represents the transpose operator.

As all receivers are situated on the Earth surface, one can consider 𝑧𝑧 1 = 𝑧𝑧 2 = 𝑧𝑧 3 = 𝑧𝑧 0 . After simplification, (5) has a form of linear equation with just 3 unknowns (𝑥𝑥 0 , 𝑦𝑦 0 , 𝑡𝑡 0 ) given by:

-2�𝐏𝐏 𝟐𝟐(𝑥𝑥,𝑦𝑦) -𝐏𝐏 𝟏𝟏(𝑥𝑥,𝑦𝑦) � ⊤ 𝐏𝐏 𝟎𝟎(𝑥𝑥,𝑦𝑦) ′ 𝑡𝑡 2 -𝑡𝑡 1 + 2𝑐𝑐 2 𝑡𝑡 0 = (𝑡𝑡 1 + 𝑡𝑡 2 )𝑐𝑐 2 - �|𝐏𝐏 𝟐𝟐(𝑥𝑥,𝑦𝑦) | 2 -|𝐏𝐏 𝟏𝟏(𝑥𝑥,𝑦𝑦) | 2 � 𝑡𝑡 2 -𝑡𝑡 1 (6)
where 𝐏𝐏 𝐢𝐢(𝑥𝑥,𝑦𝑦) represents the projection of vector 𝐏𝐏 𝐢𝐢 on the (𝑥𝑥, 𝑦𝑦) plane (i = 1, 2). In the TDoA method, time delays between different receivers are calculated considering a reference receiver. So with 3 Rx's, one can obtain two linear equations with 3 unknowns. Thus, to solve the equation for the 3 unknowns, a 4 th Rx is used. Equation ( 6) is obtained using Rx's 1-2. Similarly using Rx pairs 1-3 and 1-4, respectively, one can obtain 2 more equations. Thus, a linear system with 3 independent equations with 3 unknowns is obtained having the following form: � 𝑎𝑎 11 𝑥𝑥 0 + 𝑎𝑎 12 𝑦𝑦 0 + 𝑎𝑎 13 𝑡𝑡 0 = 𝑏𝑏 1 𝑎𝑎 21 𝑥𝑥 0 + 𝑎𝑎 22 𝑦𝑦 0 + 𝑎𝑎 23 𝑡𝑡 0 = 𝑏𝑏 2 𝑎𝑎 31 𝑥𝑥 0 + 𝑎𝑎 32 𝑦𝑦 0 + 𝑎𝑎 33 𝑡𝑡 0 = 𝑏𝑏 3 [START_REF] Wang | Passive TDOA and DOA based HF geolocation without ionosphere information[END_REF] The constants of 1 st row of ( 7) can be obtained by identification with (6) which yields:

𝑎𝑎 11 = -2(𝑥𝑥 2 -𝑥𝑥 1 ) 𝑡𝑡 2 -𝑡𝑡 1 ; 𝑎𝑎 12 = -2(𝑦𝑦 2 -𝑦𝑦 1 ) 𝑡𝑡 2 -𝑡𝑡 1 ; 𝑎𝑎 13 = 2𝑐𝑐 2 (8) 𝑏𝑏 1 = 𝑐𝑐 2 (𝑡𝑡 2 + 𝑡𝑡 1 ) - �𝑥𝑥 2 2 -𝑥𝑥 1 2 � 𝑡𝑡 2 -𝑡𝑡 1 - �𝑦𝑦 2 2 -𝑦𝑦 1 2 � 𝑡𝑡 2 -𝑡𝑡 1 (9)
Similarly from the other 2 equations obtained from Rx pairs 1-3 and 1-4, constants of other rows of ( 7) can be obtained. Finally, the transmitter coordinates (𝑥𝑥 0 , 𝑦𝑦 0 ) and the transmit time 𝑡𝑡 0 can be obtained by doing a simple matrix inversion operation. The value of 𝑧𝑧 0 ′ can be obtained by substituting (𝑥𝑥 0 , 𝑦𝑦 0 ) in [START_REF] Torrieri | Statistical theory of passive location systems[END_REF]. Finally the equivalent height ℎ 𝑒𝑒𝑒𝑒𝑒𝑒 can be obtained by resolving [START_REF] Pagani | Time domain HF geolocation: Experimental measurements and preliminary results[END_REF].

Note: Equation ( 7) can be realized in the form of (𝑥𝑥 0 , 𝑦𝑦 0 , 𝑡𝑡 0 ) ⊤ = 𝐀𝐀 -1 𝐁𝐁 to obtain the transmitter coordinates. In the case of more receivers, the matrix inversion is done using the pseudo inverse of matrix A.

IV. SIMULATION RESULTS

A. Quasi-Parabolic Ionosphere with a Fixed Height Fig. 3 represents geolocation error in locating Tx on a square grid with a length of 4000 km using the least number of 4 receivers required by the TDoA method. The 4 receivers are placed in a manner that they approximately form a square with the center at the origin of the square grid. The artifact seen along both the axes of the grid in Fig. 3 is due to the approximate symmetry formed due to the manner in which the receivers are placed. The length of the sides of the square formed by the receiver is approximately 420 km. These four receivers are illustrated by green "+" symbols in Fig. 3. The group path lengths between different receivers and all possible transmitters on the grid are evaluated using the analytic expression of the QP model. The values used to define the QP model in all simulations are: F = 0.8, hm = 300 km and y m = 100 km. The coordinates of the transmitter can be found by the intersection of three hyperboloids using the geolocation algorithm of TDoA explained in Section III. The height of maximum ionization (ℎ 𝑚𝑚 ) of the QP ionosphere is assumed to be fixed, thus leading to a uniform ionosphere. The geolocation algorithm being dependent on the reference receiver, the Tx coordinates are evaluated considering each receiver as the reference receiver. The final Tx location is obtained by computing the mean of the Tx coordinates obtained from each reference receiver. Finally, the geolocation error is evaluated by calculating the difference between the true location and the obtained final location. All the geolocation error greater than 100 km is represented by black color.

The observed geolocation error is very large in the case when the Tx site is located at a distance less than around 600 km from a Rx site. This geolocation error is higher due to large variation of the equivalent height between a Tx and different Rx configurations when they are placed close to each other. However, as a Tx site is further away from the origin, the geolocation error decreases. The Tx located at a distance range of around 1200 to 2000 km from the origin of the square grid could be located with an error of about 20-50 km using four receivers.

Fig. 4 represents geolocation error calculated for each Tx location using five receivers with the TDoA method. The 4 receivers are placed in a way that they approximately form a square with the length of side approximately equal to 710 km and the fifth receiver is placed at the centre of the square grid. The geolocation error for each Tx location, assuming refraction from a uniform QP ionosphere is calculated in a similar manner as done in the four receiver case. As seen in Fig. 4, the maximum geolocation error is around 50 km very close to the receivers. With the use of an additional receiver and placing it in an optimum position, the geolocation accuracy is improved.

Fig. 5 represents the CDF of geolocation error in locating transmitters on a square grid with a length of 4000 km using different configurations of the number of receivers which are placed approximately uniformly over the grid. It is also assumed that the refraction of radio signals is from a QP ionosphere with fixed height. The numbers in the legend box of Fig. 5 corresponds to the number of receivers. From Fig. 5, it is also evident that 80% of the geolocation error is less than 95 km using four receivers whereas 80% of the geolocation error is less than 30 km using five receivers. With five receivers, the median of the geolocation error is reduced by almost 85% compared to the median of the geolocation error with four receivers. This gain in terms of reduced geolocation error is quite high at the expense of just one additional receiver. Furthermore, with increasing the number of receivers the geolocation error keeps decreasing but the improvement is less significant. As seen from Fig. 5, there is only a slight improvement in the geolocation error when the number of receivers is increased from 13 (yellow curve) to 17 (cyan curve). From Fig. 5, it is also evident that 93% of the geolocation error is less than 10 km and the maximum geolocation error is just 25 km using 17 Rx's. Note that the Tx coordinates is evaluated using only the receivers for which the transmitter lies within their radio horizon.

B. Quasi-Parabolic Ionosphere with a Variable Height

In order to obtain the location of the Tx using more realistic conditions of the ionosphere, the height of maximum ionization ℎ 𝑚𝑚 is varied as a Gaussian distribution with a mean of 300 km and variance of 20 km all across the square grid of length 4000 km. Fig. 6 represents the CDF of geolocation error in locating transmitters on a square grid with a length of 4000 km with different configurations of the number of receivers and the height ℎ 𝑚𝑚 being varied as a Gaussian distribution. The numbers in the legend box of Fig. 6 corresponds to the number of receivers. As expected, the geolocation error is larger in this case for a similar number of receivers when compared to refraction from a uniform QP ionosphere. From Fig. 6, it is seen that around 20% of the geolocation error are less than 100 km using four receivers whereas around 75% of the geolocation error are less than 100 km using five receivers. The geolocation accuracy improves considerably when the number of receivers is increased from 5 to 9 and then, from 9 to 13 as shown in Fig. 6.The geolocation error keeps decreasing with increasing number of receivers.

V. CONCLUSIONS

This paper presents, HF geolocation using TDoA with variable number of receivers. These techniques make use of single antenna receivers which are a cost effective solution and can be easily deployed. The group path between a Tx-Rx pair for oblique propagation is obtained using the QP model. Furthermore, the mathematical approach for geolocation is also explained for the TDoA case. Simulation studies show that the geolocation error is significantly reduced by increasing the number of receivers. 95% of the geolocation error is approximately less than 10 km using 17 Rx's for a uniform QP ionosphere; 95% of the geolocation error is less than 40 km using 17 Rx's for a QP ionosphere having a variable height. The position of the receivers also plays an important role in the geolocation of the Tx. Several aspects related to the position of the receivers will be investigated in our future research. An experiment will be setup in the near future with four or five receivers to confirm these simulation results.

Fig. 1 .

 1 Fig. 1. Plot of ground distance using QP model at different values of F.

Fig. 2 .

 2 Fig. 2. Geometry of the HF localization principle considering uniform single layer ionosphere.

Fig. 3 .

 3 Fig. 3. Geolocation error using TDoA with 4 Rx's for a QP ionosphere.

Fig. 4 .

 4 Fig. 4. Geolocation error using TDoA with 5 Rx's for a QP ionosphere.

Fig. 5 .

 5 Fig. 5. CDF of geolocation error for different numbers of Rxs assuming a QP ionosphere with fixed height.

Fig. 6 .

 6 Fig.6. CDF of geolocation error for different numbers of Rxs assuming a QP ionosphere with variable height following a Gaussian distribution.

ACKNOWLEDGMENT

Part of this study has been funded by the Carnot Institute and the Brittany Region of France.


Justin Berli
email: justinberli@gmail.com

Mattia Bunel
email: mattiabunel@hotmail.fr

César Ducruet

Sea-land interdependence in the global maritime network: the case of Australian port cities

Keywords: accessibility, Australia, centrality, complex networks, intermodalism, maritime transport, multiplex graph, port cities, shipping flows

This article tackles the longstanding issue of intermodality head on. From a geomatics perspective, we model both maritime and road networks connecting port and non-port cities taking into account crucial features such as physical geography, shortest paths, and transport costs. This creates the opportunity to study a hybrid networkboth planar and non-planar, and the centrality/accessibility of cities in this bilayered network. Based on the case of Australia, main results convey new empirical findings on how port and urban hierarchies correlate with single-layered and bi-layered connectivity. We discuss main results in the light of network science, spatial science, and transport studies.

Introduction: background and objectives

Measuring and analyzing intersections between networks is a major challenge for both researchers and practitioners. In transport geography and economics, there is a consensus that the efficient combination of sea and land accessibility is necessary -if not vitalto port activity and to the expansion of wider supply chain systems. Yet, this perspective still lacks rigorous empirical validation today, taking into account the topological structure of both land and maritime networks. Early geographers long focused on ports and port cities at the local level of the port-city interface from a morphological perspective [START_REF] Rimmer | Maritime Networks. Spatial Structures and Time Dynamics[END_REF]Robinson, 2015a). During many decades, numerous port studies (see [START_REF] Ng | Port geography at the crossroads with human geography: Between flows and spaces[END_REF] for a synthesis) left behind continental and maritime connectivity issues despite the existence of conceptual frameworks such as the "port triptych" back in the late 1970s [START_REF] Slack | Epilogue[END_REF], in which the port is a node connecting hinterland (inland market area) and foreland (overseas market area). Considering these two objects as one single, intermodal and interconnected, network, however, is subject to a number of issues deserving discussion across various research fields. One of them is the urban perspective on global networks [START_REF] Taylor | World City Network. A global Urban Analysis[END_REF], which mainly focuses on multiple immaterial linkages created by multinational firms, with less importance given to transport networks and especially to ports [START_REF] Jacobs | Integrating world cities into production networks: The case of port cities[END_REF]. In geography and elsewhere, mainly local/national planar and technical networks have long been at center stage [START_REF] Ducruet | Network science and spatial science: Review and outcomes of a complex relationship[END_REF] until the extension of such analyses to non-planar transport and global networks. Yet, transdisciplinary urban studies keep on insisting about the need to understand the metabolism [START_REF] Kennedy | Energy and material flows of megacities[END_REF] and the "pace of life" of cities [START_REF] Bettencourt | Growth, innovation, scaling, and the pace of life in cities[END_REF], understanding that efficient connectivity is favorable to urban growth (see also [START_REF] Neal | The urban metabolism of airline passengers: Scaling and sustainability[END_REF].

Inter-network externalities

The analysis of combined networks witnessed a growing emphasis from researchers during the last decade or so. Changes in topological structure inferred by the interconnection of different networks (and/or different layers within a single network) occupy the bulk of related research (see Vespignani, 2010) on both nodes and links. Although analyses of single networks remain dominant in the literature, it is possible to find several examples of works looking at least at two types of linkages. Throughout the latter perspective, most empirical studies looked at networks or layers of the same nature: planar networks, such as electrical transmission network and Internet network backbone [START_REF] Rosato | Modelling interdependent infrastructures using interacting dynamical models[END_REF] or non-planar networks, such as airline networks and corporate networks [START_REF] Liu | Exploring the co-evolution of the geographies of air transport aviation and corporate networks[END_REF], to name but a few (see a recent review by [START_REF] Ducruet | Multilayer dynamics of complex spatial networks: the case of global maritime flows (1977-2008)[END_REF].

The present research wishes to investigate the interconnection between two networks of different structure, namely planar and non-planar networks. Planar networks can be defined as two-dimensional graphs where links cannot cross without creating a node at the place of the intersection. This is the case of more land-based transport networks but also of so-called "technical networks" in general from electric power grids to the Internet. In geography, one of the earliest examples of the kind was provided by [START_REF] Cattan | Barrier effects: The case of air and rail flows[END_REF] in her study of barrier effects affecting air and rail flows between European cities (see also [START_REF] Van Geenhuizen | Interconnectivity of transport networks: A conceptual and empirical exploration[END_REF]. Later on, [START_REF] Berroir | Les systèmes urbains français[END_REF] revealed the emergence of subnetworks among French urban areas based on multiple interactions such as air flights, commuting flows, scientific collaborations, high-speed rail linkages, and firms' networks. This particular approach of the "French school" (P.A.R.I.S. team) takes its roots, at least in part, in earlier works of the so-called "New Geography" on the Swedish side (see a review by [START_REF] Peris | Penser les villes en réseaux : une analyse des théories sur les liens interurbains[END_REF], with the pioneer works of the geographer G. Törnqvist, for instance, who proposed already to study both material and immaterial flows linking cities. More recently, generative network models were applied to a composite transport network defined as a "surrogate measure of three individual transport networks: road, rail, and air transport" among Southeast Asian cities (Dai et al., 2016). Recent attention has been paid, however, to the importance of road transport for supply chain performance using a multiplex network framework [START_REF] Viljoen | The road most travelled: The impact of urban road infrastructure on supply chain network vulnerability[END_REF]. Such studies have in common, however, to neglect the importance of shipping linkages, depending on the context of the study area, and to focus on a composite network. The latter perspective does not allow grasping the respective role of the individual networks or layers in the overall centrality of nodes or cities. The main purpose of our paper is thus to untangle the effects of each network or layer and to better understand the specific role of shipping linkages in the connectivity of cities.

Intermodal transport and port cities

The fact that ports and port cities benefit from being well connected to both sea and land networks is common thought in transport studies and port geography, referring back to the port triptych model proposed by the French geographer André [START_REF] Vigarié | Ports de commerce et vie littorale[END_REF] in which ports connect foreland and hinterland. Recurrent arguments include the fact that sea-land transshipment may generate economies of scale while increasing the number of destinations and shipping frequency, thereby opening new markets to ports beyond their traditional hinterlands (van Klink and van den [START_REF] Van Klink | Gateways and intermodalism[END_REF]. While the influence of physical infrastructure interconnection cannot be ignored [START_REF] Bottasso | Port infrastructures and trade: Empirical evidence from Brazil[END_REF], most studies focus more on the role of vertically integrated actors in establishing intermodal services [START_REF] Franc | Analyzing hinterland service integration by shipping lines and terminal operators in the Hamburg-Le Havre range[END_REF][START_REF] Notteboom | The relationship between port choice and terminal involvement of alliance members in container shipping[END_REF] and so-called dry ports (inland) to relieve seaports from congestion and enhance distribution systems, among other goals [START_REF] Roso | The dry port concept: Connecting container seaports with the hinterland[END_REF]. Hinterland delineation, a longstanding issue in transport studies, recently improved but tended to ignore shipping linkages per se [START_REF] Halim | A strategic model of port-hinterland freight distribution networks[END_REF][START_REF] Tiller | Spatial patterns of landside trade impedance in containerized South American exports[END_REF][START_REF] Moura | Delimiting the scope of the hinterland of ports: Proposal and case study[END_REF][START_REF] Jung | Community structures in networks of disaggregated cargo flows to maritime ports[END_REF]. Most related empirical investigation of sea-land interaction remains qualitative to date (see [START_REF] Torbianelli | When the road controls the sea: a case study of Ro-Ro transport in the Mediterranean[END_REF]. Such an approach is more developed in airline network studies such as [START_REF] Choi | Comparing world city networks: A network analysis of Internet backbone and air transport intercity linkages[END_REF] crossing airline networks and the Internet backbine, [START_REF] Matisziw | Evaluating locational accessibility to the US air transportation system[END_REF] looking at air transport flows in relation to land access, and [START_REF] Derudder | The connectivity of South Asian cities in infrastructure networks[END_REF] combining airline flows with planar networks such as bus schedules, railways and roads in Southeast Asia. In fact, across a sample of 728 scientific articles published about ports in geography journals and other disciplines between 1950 and 2012, [START_REF] Ng | Port geography at the crossroads with human geography: Between flows and spaces[END_REF] found that those dealing with supply chain linkages, intermodal transportation, and networks rapidly increased in number, although very few of them actually offered an empirical, quantitative analysis of combined sea-land networks connecting ports.

An early breakthrough in that direction by [START_REF] Chapelon | L'accessibilité, marqueur des inégalités de rayonnement des villes portuaires en Europe[END_REF] looked at the inland accessibility of European seaports, simulating their respective population and Gross Domestic Product potentials based on trucking time, road quality, driving regulations, tolls, etc. One originality of his study was to include a handful of ferry linkages to complement the sole road network, but the European maritime network remained too simplified to account for a realistic account of shipping connectivity. A radically different approach came out with the work of [START_REF] Nelson | Travel time to major cities: A global map of accessibility[END_REF] providing a composite index of accessibility for the world's cities. No less than twelve layers compose this index using raster methods, including AIS data for shipping flows (Automated Identification System), road and railway networks, elevation, land use, etc. [START_REF] Halim | The impact of the emergence of direct shipping lines on port flows[END_REF] used Spatial Computable General Equilibrium (SCGE) modelling to explore the impact of opening new shipping routes on port activity, including a simplified road network connecting ports with hinterlands. More recently, [START_REF] Guerrero | Foreland mix and inland accessibility of European NUTS-3 regions[END_REF] provided a new analysis of sea-land connectivity using subnational regions as units of analysis, AIS data, and Euclidian distances between ports and inland markets, while Shen (2017) combined trade and customs data with road transport networks to analyze door-to-door shipping flows between U.S. and Chinese cities through GIS methods. In a similar vein than [START_REF] Chapelon | L'accessibilité, marqueur des inégalités de rayonnement des villes portuaires en Europe[END_REF], [START_REF] Wang | Spatial and economic effects of the Bohai Strait Cross-Sea Channel on the transportation accessibility in China[END_REF] explored the effects of opening a new coastal shipping artery across the Bohai Rim on local economic development and sea-land accessibility. Other interconnection studies combined shipping linkages with other networks such as airlines (Parshani et al., 2010) and trade flows [START_REF] Catalayud | Connectivity to international markets: A multi-layered network approach[END_REF].

What remains lacking is the modeling of the shipping network and the analysis of its interconnection with landside networks, especially in an urban context. Most of the time, network nodes are transport terminal and/or ports thereby ignoring their urban context and dimension (see [START_REF] Qu | Sustainability SI: Multimode multicommodity network design model for intermodal freight transportation with transfer and emission costs[END_REF] on an intermodal study focusing on greenhouse gas emissions). Although the contemporary literature dominantly insisted on the dereliction of port-city linkages, several studies argued in favor of maintained positive externalities offered by cities to port activities [START_REF] Hall | Why are maritime ports (still) urban, and why should policy makers care?[END_REF], demonstrating that the global maritime network had remained tied to cities and, in particular, extended city-regions crossing vessel traffic and demographic size as the main indicators over the last 120 years [START_REF] Ducruet | Maritime networks as systems of cities: The long-term interdependencies between global shipping flows and urban development (1890-2010[END_REF]. This is the goal of the present research to tackle such lacunae and to go one step further by combining the road network with the maritime network to study urban functions and hierarchies.

Main objectives and application to the Australian case

Despite all the major advances made by scholarly research on sea-land connectivity, a number of issues remain unresolved and/or unexplored. Firstly, it is rather common that the urban dimension of network nodes is largely ignored, the terminal or port (or the subnational region) being the main unit of analysis. Thus, it remains unclear whether sea-land connectivity influences, or is determined by, urban activity (e.g. population) or port activity (e.g. throughput). This paper thus adopts a new perspective, considering how urban hierarchy and port hierarchy correlate with mono-layered and bi-layered centrality. This has implications for both network studies, transport studies, and urban studies, in terms of specialization and multiplexity. Second, the concrete spatial distribution of the shipping network remains ill defined, i.e. without taking into account its geographical and operational features. This article proposes a new approach through the modeling of a maritime grid capable of mapping and measuring shipping connectivity in a more precise manner, especially in relation to landside transport networks.

The choice of Australia is motivated by several aspects. First, it is an island country, which avoids hinterland competition like in Europe or USA among ports serving the same "contestable" markets, as Australian port cities all have their own port and there is limited traffic shifts among them. Population is concentrated along the coastline in those cities mainly, so that the heartland of the country is relatively empty. Second, Australia as a whole underwent steady vessel traffic growth since the late 1970s (Figure 1), with a rising (albeit relatively minor) world share, thus motivating the study of its maritime connectivity. Currently, Australia is the world's 23 rd largest export economy, although its imports surpass its exports in terms of trade value, exporting mainly solid bulks (iron ore, coal briquettes, gold, wheat, and crude petroleum) and importing liquid bulk (refined and crude petroleum) as well as finished goods (cars, computers, and packaged medicaments (Observatory of Economic Complexity, 2018). It is now famous for having become the world's largest exporter of iron ore and the world's leading coal exporter, partly fostered by Chinese demand [START_REF] Wang | Transport corridors and regional balance in China: The case of coal trade and logistics[END_REF].

[Insert Figure 1 about here] Australia's international trade thus occurs mainly by sea up to 99% of its total, and is highly diversified serving a developed economy and serving both specific industries and a population of about 24 million inhabitants concentrated in large coastal (port) cities. Heavy industries need considerable port infrastructures close to the mines located mainly in the Northwest of the country for their exports, while main urban centers locate more in the Southern part, from Perth in the West to Sydney in the East and Brisbane more in the Northeast. Those latter locations are mainly responsible for the trade and import of manufactured goods and for passenger traffic as well as oil imports to feed the urban system with energy. The Australian port system, and especially bulk ports, thus necessitates efficiency intermodal facilities to connect production and consumption centers to global markets (see [START_REF] Robinson | Regulating efficiency into port-oriented chain systems: export coal through the Dalrymple Bay Terminal, Australia[END_REF]2015b). As such, this port system is marked by a high degree of specialization (Figure 2) whereby a majority of ports handle solid bulks and the rest is specialized in one main traffic category, such as container ports being also large cities: Melbourne, Adelaide, Port Botany (Sydney), Brisbane, and Fremantle (Perth). Except from liquid bulk that is still handled by numerous ports, a minority is specialized in one main traffic type, such as passengers including cruise and ferries (Hayman Island, Broome, Exmouth, Eden, and Burnie), food products including fishing and drinks (Darwin), and general cargo (Port Alma and Port Lincoln). Although the figures might have changed since the study year 2008, there is a high stability in the type of products handled by either large cities or specialized bulk ports. This look at the Australian port system in relation with trade and settlements complements recent works focusing on one particular port city such as Sydney, looking at traffic specialization and fluctuation overtime [START_REF] Paflioti | Box dynamics: A sectoral approach to analyse containerized port throughput interdependencies[END_REF]. Other earlier studies investigated the strong influence of cities (rather than ports) in the attraction of the maritime service industry in Australia [START_REF] O'connor | Australian ports, metropolitan areas and trade-related services[END_REF]) and later at the world scale [START_REF] Jacobs | The location and global network structure of maritime advanced producer services[END_REF], long after pioneering studies of Australian ports where urban and hinterland elements where mentioned [START_REF] Rimmer | The search for spatial regularities in the development of Australian seaports 1861-1961/2[END_REF] to explain at least partly port dynamics in reference to classic models of port system evolution. Earlier works looked, for instance, at the rank-size distribution of coastwise shipping flows among Victoria state ports [START_REF] Britton | Coastwise external relations of the ports of Victoria[END_REF].

[Insert Figure 2 about here] The remainder of the paper are as follows. Next section 2 presents the methodology for modeling and combining shipping and road networks, focusing on the particular case of Australia as a testing ground. Section 3 proposes computes and compares classic centrality/accessibility algorithms on individual and combined networks, while section 4 identifies the degree to which such measures correlate with port and urban hierarchies in Australia. The final section serves as a discussion and conclusion about the contribution of the findings to existing research and opens pathways for further research in the field of urban network studies.

Modeling the sea-land transport network

2.1 The worldwide maritime grid

Preliminary specifications

Visualizing maritime flows is a longstanding issue in geography and other sciences, although the digital revolution of the 2000s resulted in a variety of solutions and experiments, most of them being accessible online [START_REF] Bunel | Geovisualizing the sail-to-steam transition through vessel movement data[END_REF]). Yet, most research in this domain remains dedicated to radar or satellite data that are more precise than vessel movements, which only indicate ship movements between ports of call. Our proposal in this paper is to deliver a new and robust architecture notwithstanding the necessity to overcome numerous development issues, based on the steps described in Figure 3.

[Insert Figure 3 about here] A maritime grid thus provides an approximation of paths taken by ship movements while respecting the geographic constraints of coastlines. This allows defining how to segment and partition the continuous ocean space where ships circulate by using a virtual grid. Our experiments consisted in creating a first maritime grid of constant mesh [START_REF] Bunel | Geovisualizing the sail-to-steam transition through vessel movement data[END_REF]. But this choice quickly turned out to be irrelevant. Depending on the mesh chosen, the grid was either too inaccurate or too dense to be easily exploitable. In front of this difficulty of finding a balance between these two opposite points, we develop a new methodology, based on the observation that the accuracy of the grid did not have to be constant (Figure 4). Our model rests on a regular meshing worldwide composed of eight squares of 90°side. From this starting element, the meshing was refined with an iterative process. We subdivided each square intersecting a continent into four same sized squares, except for those fully included in a continent or an ocean. Each iterative step multiplied the number of squares, their area being potentially divided by four, resulting in a denser meshing near coastlines. The resulting mesh is strongly dependent of the coastal line and, by expanding, of the data used. In our case, we used data from the Natural Earth project (with a 50m scale). In this case, seven iterations, providing the best ratio between complexity and accuracy, resulted in a final meshing composed of 23,000 squares with areas extending between 6 x 10³ km² and 6 x 10⁶ km² (Figure 4a). This grid is the basis of ulterior calculations. A trajectory is defined by a list of adjacent squares, computing the shortest paths (i.e. using pg_routing, a PostGIS extension providing geospatial routing functionality to the database) through linking the squares' respective centroids (Figure 4c), using the Moore neighborhood (two squares are neighbors if they share at least one vertex) (Figure 4e).

[Insert Figure 4 about here]

The created grid allowed a fluid mapping of trajectories on seas and oceans while including rivers and lakes. It was completed by new links based on Natural Earth physical data. Once rivers and lakes links were connected to the grid, the final procedure was the creation of links between the grid and the actual position of ports. The first attempt was the linking of each port to the closest node of the grid. However, this solution was not exploitable: the created links were intersecting continents. Closest ports were thus grouped into clusters, and those cluster centroids were linked to the grid. Finally, the links were weighted according to the distance between their two extremities. Connected components deriving from the aforementioned steps were then linked with the nearest and largest connected component to avoid disconnecting the network, such as in the cases of closed seas (e.g. Black Sea), straits (e.g. Gibraltar), and canals (e.g. Suez, Panama). Using this set of weighted links, we are able to calculate an approximation of the shortest sea road between two ports (and more generally between any pair of points). Although necessarily imperfect, this approximation remains effective and allows us to reconstruct coherent maritime trajectories with those observed with recent data, such as satellite images.

Maritime flow data

The shipping data used in this study was obtained from Lloyd's List, the world's main maritime insurance company and marine intelligence. The database consists in the aggregation of daily vessel movements occurring during the complete months of January, March, July, August, and September 2008. The trajectory of each vessel was computed at the level of links and of port nodes in the form of a global, weighted origin-destination matrix including direct and indirect voyages among ports of the world, based on 920,014 movements made by 46,370 vessels of all types connecting 3,981 ports. With reference to Hu and Zhu (2009), direct voyages correspond to the space-L topology where nodes are connected by inter-port flows without calling at an intermediary node, while indirect voyages correspond to the space-P topology where we add all inter-port links having calls at least at one intermediary port. The inclusion of indirect linkages allows considering longer-distance linkages; otherwise a network made of direct links only would not connect Australia with ports outside the Asia-Pacific area. In other words, the degree centrality of Australia is limited to direct neighbors in space-L but it expands to all connected ports in space-P, with a depth varying according to the length of voyages. For instance, a vessel voyaging from Melbourne to Singapore and then Dubai would not connect Melbourne and Dubai in space-L (only two links) contrary to space-P that would also container this third link.

Based on the interconnected sea-land network, we propose measuring the centrality of Australian port cities based on classic measures borrowed from graph theory and complex networks, namely Shimbel Index (i.e. number of links between each node and its farthest node) transformed into its inverse value, closeness centrality, and betweenness centrality (i.e. number of occurrences of nodes on shortest paths). Vessel movement data from the Lloyd's List were stored in a PostGreSQL/PostGIS database. This database is requested by Django, a Python web framework. The Python choice as the principal language of development on the server side (i.e. for the calculation done on the database side) allows the use of numerous libraries, specialized in scientific calculation (Numpy), or graph treatment (NetworkX). By opting for this choice, the user is able to visualize different types of statistical or graph-theoretical indicators on chosen data. On the client side (i.e. on the web browser), the results are mapped with OpenLayers and Cesium javascriptJavaScript libraries, completed by graphs done with D3.js.

The road transport network

Preliminary specifications

The road network model is currently being developed for OpenStreetMap (OSM) data as it is the most precise, up-to-date database available to conduct studies on the 2000/10 period. The precision, the completeness of the main road system (Hakley, 2010;Wang et al., 2013), and the collaborative nature of this dataset provide the best solution to create a coherent network. Given the level of detail provided by this data source, some lightering operations must be deployed to generate a network for the entire world in a reasonable amount of time and without reaching computational limitations. Throughout this section, all presented features and processes are made and run automatically using a Python script combined with PostGIS queries. The objective of the model is to create, for the entire world, a routable network consisting of road sections, connecting portsthose extracted from the Lloyd's Shipping Index and other sources (cf. previous section) and main cities. Australia has been chosen mainly because it presents an enclosed road network and dense cities as well as barren, sparsely inhabited lands. This choice is, in a way, purely arbitrary and it means that the model will have to be enhanced to fit the diversity of other territories.

Being able to compute the entire world is an objective that induce the use of a partitioning phase. Indeed, some spatial queries can be highly time-consuming andas of the huge number of road sectionsbe mobilizing a significant amount of physical resources. Furthermore, running a unique query for the entire world can provoke server-sided memory issue and would need days of processing; any raised error would then require the query to be resumed, thus, significant amount of time to be lost. The most suitable partitioning solution turned out to be the quadtree, which allowed us to create tiles based on the density of cities and ports. This method recursively divides the world into four areas of identical extent and repeats itself until the number of entities per tile is below a user-defined value. The quadtree is useful to create homogeneous tiles in terms of ports and cities, and, consequently, of the associated road network. The main issue raised by the use of any partitioning methods is the loss of connectivity as each area is independently processed. Road sections crossing the tiles boundaries must be taken into consideration and handled specifically to ensure the creation of a consistent, connected network. This phase increases the amount of processing time but is the inevitable counterpart of the use of spatial partitions.

For now, the model is divided into four distinct stages which are essential given the fact that some contextual operations need to rely on data already computed in the surrounding tiles. Each tile are processed four times, and during each stage, a different operation is performed: first, littoral areas are derived based on the primary road system ; second, a contextual extraction is applied on littoral road sections to ensure that each ports are connected to the primary road system ; third, a routable network is created with informations such as source node, target node and the cost ; fourth, this network is then simplified to keep only ways used by Dijkstra's shortest paths algorithm. Every stage will be explain with more precision in the present section.

Data extraction

Before constructing the actual network, some density-reducing operations need to be performed to get rid of unwanted road sections (Figure 5). This will fasten the creation of the network mainly because less roads will be taken into consideration during this operation which imply the use of lots of spatialrelated functions. First, the creation of littoral areas based on the primary road network will allow the secondary road system to be extracted. Then will be kept only littoral road sections around ports depending of an internally defined hierarchy. The primary road system is defined by the use of the highway=motorway and highway=trunk from the OSM dataset as it gives a consistent overview of the main road network connecting important cities [START_REF] Schmid | A Method for High-Level Street Network Extraction of OpenStreetMap Data in OpenScienceMap[END_REF]. This constitutes a crucial stage in the model as it will significantly reduce the amount of data to process by getting rid of unnecessary road sections. The main idea lies in the creation of littoral areas based on the primary road system. Those areas are not actual littorals, they are internally defined, as it only depends on roads and does not take into consideration geomorphological features, climate, vegetation, etc. First, the primary road sections of the entire world are stored inside a temporary table. Then, for each tile, every road located inside the surrounding tiles are transformed into polygons, which are kept only if they overlap the sea, and finally are merged into a unique entity if they have any border in common. Those polygons are formed using the road system and the tile boundary combined. This method allows us to define an area in which ports are likely to be located. It is then possible to extract the secondary road sections (tags highway=primary, secondary, tertiary, residential, service, road, track, unclassified from OSM) inside those littoral areas, consequently ensuring that every ports will be connected to the road network and reducing the amount of data to be processed. Processing the surrounding tiles seems to be the only solution to obtain seamless littoral areas throughout the world, even though the duration of the query is increased.

Following the extraction of littoral secondary roads, another extraction operation needs to be performed to once again get rid of unnecessary data. This method can be qualified as contextual as it is based on ports, the main object of the study, and relies on a model-intrinsic hierarchy that defines a distance in which road sections will be kept. First, a rank is assigned to each ports, depending on the most important road located nearby. That rank will determine the importance of the expected road network that connects the port to the continental primary road network. Then, the Euclidean distance between the port and the primary road system is calculated, multiplied by 1.5 and used as a buffer area in which only roads of equivalent or higher rank are extracted. This method implies that, for instance, a tertiary port (i.e. located nearby a tertiary road) is connected to the primary road network only by tertiary, secondary and primary road sections.

[Insert Figure 5 about here] Theoretically, this approach is consistent regarding the overall network structure but, given the nature of the data and some particular features specific to a given territory, some issues are notable afterwards and justify handmade corrections. Sometimes indeed, a low-rank road section links two of higher, therefore destroying the network integrity as each portclassified according to the closest roadis linked to the main road system considering surrounding higher rank sections in order to reduce data density. Currently, those problematic situations are still sparse and can be easily handled by the operator. Another issue is a direct consequence of partitioning the data using the quadtree (cf. 1.2.) and fixing it was accomplished by considering the surrounding tiles when extracting road sections linking a port to the rest of the network. Then, road sections from adjacent tiles were inserted into the table if they ensured the connectivity of the ports to the primary road system; naturally, after checking if those road sections hadn't already been inserted.

Only continental primary and littoral secondary road sections have henceforth been extracted to lighten the highly dense OSM dataset. The creation of a routable network based on those roads requires the use of queries that inherently demand a certain amount of resources, specifically those deploying spatial relationships between geographical entities (e.g. if they intersect, overlaps, touches, etc.). The goal of this phase is to obtain a table containing every extracted road with specific routing attributes (source node, target node, cost), calculate shortest paths between each cities and ports and finally, keep only road sections that are used by the algorithm.

Base network computation

The computation of the base network uses the combination of the continental primary and littoral secondary road sections. On that basis, start and end nodes of each sections are extracted and used as a point for the roads to be split, thus allowing to retrieve every road junction in the network. Then, the entire network is dissolved without bridges (i.e. lines that crosses each other) to obtain a multipolygon which is then divided into sections defined by intersection points and converted into simple polygons. A speed value is mandatory to calculate costs; the problem is that the previous operation got rid of every attributes. The idea is to perform a spatial join based on a containing clause; if the new road section contains the former ones, then the join is made (Figure 6). Finally, costs are calculated using a simple formula: cost (km/km/h) = length (km) / max speed (km/h).

[Insert Figure 6 about here] Another operation is needed to add information and transform this dataset into a routable network. Indeed, the routing algorithm requires the source and the target node of every section to be able to navigate through the network system. The nodes are extracted simply by calculating distinct start and end points of every section and by assigning them a unique identifier. Then a spatial join is performed based on the intersection between the start and end points of each road section and the newly-created nodes. In the end, the network obtained is routable, bridges are still present and every section is weighted according to its speed limit and its length.

Having a routable network now would be enough for the implementation. The problem is that the density-reduction operations, by their structure, extracted and kept residual road sections that are superfluous, as they are not used to link ports to the continental primary road system. A final stage is then essential to eliminate them and guarantee a fully lighten network: calculating the shortest paths between each main nodes of the network (ports and cities) and keep only necessarily followed road sections. To achieve this goal, the PostGIS implementation of Dijkstra's algorithm is deployed and are extracted only roads sections that are used by the computed shortest paths. Keeping the connectivity between tiles can be ensured by considering all surrounding tiles during the process. Secondly, the road sections that have been extracted are clipped according to the tile boundaries so that no overlapping features will appear in the final table. The network is re-calculated the same way as in 3.1. to dump intermediary nodes, artifacts of the previous network computation. The only unnecessary nodes kept are those located onto the tiles boundaries. Finally, source and target attributes of the road sections are updated according to the new network nodes (Figure 7).

[Insert Figure 7 about here]

The result for Australia is realistic as in addition to the port traffic hierarchy, our methodology provides a view of the related maritime network and how those ports are at the same time, connected to the modeled road network (Figure 8). One interesting fact is that the main maritime routes connecting Australia with its neighborhood (e.g. Malacca Strait) are concentrated along the east coast where main urban centers are located, with routes conducting to Northeast Asia and Malacca / Indian Ocean, while the largest bulk ports in the North remain connected through lower traffic links despite a comparable port traffic than those cities. It also reflects the trade imbalance of Australia as imports surpass exports in value but also in volume.

[Insert Figure 8 about here]

Intermodal accessibility: application to Australian port cities

Calculating intermodal centrality

To highlight the specificities of Australian ports' accessibility, the combination of both maritime and road networks has been deployed to calculate betweenness and closeness centralities at different levels. A simple process to create an intermodal network consisted in combining terrestrial roads and maritime potential routes. Several distinct cases have been used to calculate indexes, thus allowing a better overview of the ports' accessibility. They all directly concern the maritime grid, depending on how it is handled, and represent different abstraction levels. The first case considers the complete maritime grid, which represents potential routes that vessels can use; the second case considers an Origin-Destination (O-D) matrix for 2008, which ignores spatial constraints (distance, continents, etc.) to keep only a boolean value between ports (presence or absence of a link); the third is a blend of the two previous methods: it affects shipping flows in 2008 to on the complete maritime grid, thereby adding spatial constraints. On the other hand, the road network is always taken in its integrity. On that basis, for every of those three cases, betweenness and closeness centrality indexes have been calculated considering both the combined sea-land network, and each network separately. Running accessibility indexes calculation on a spatial network requires a weight value for every edge, as we compute a shortest paths algorithm throughout the graph. The main issue is that the road network's edges are weighed by their cost, calculated according to their length and the speed limits (cf. previous section). Without proper information about vessel speed, it was impossible to apply similar methods to weight the maritime grid. Therefore, only indexes calculated on the sole road network use the speed-derived costs while those calculated on the maritime grid or on the combined sea-land network use only topological distance. This approach remains consistent, as the calculation of indexes applies to a homogeneously weighted graph. The actual calculation of betweenness and closeness centrality indexes were performed using the free programming language R and its Igraph network analysis' library collection. For now, only indexes concerning ports are interesting as the road network cannot be considered as a proper graphit does not contain nodes (cities and ports) directly linked by edges but a sum of small routable sections. In some cases, a shortest path may run through a city's ring highway and not by the actual city node, thus inducing wrong centrality values. Ongoing efforts are made to remedy such a "noise", by removing intra-urban road sections to retain only intercity road segments.

Centrality scores and land-sea interactions

Landside issues in sea-land operations of Australian port cities is a recurrent concern for scholars and related supply chain actors [START_REF] Robinson | Port-oriented landside logistics in Australian ports: A strategic framework[END_REF]. Results obtained from the computation of accessibility or centrality measures on single or combined networks provide interesting evidences about their respective roles for Australian port cities (Table 1). One main result is that closeness centrality is highly correlated for both island and mainland ports, whatever the supporting network, but mainly when comparing maritime centrality and combined centrality. The exception is for mainland ports using the maritime grid as the base network as correlation is near zero. Perhaps, the high density of the grid compared with the other networks had over-or underestimated Australian port's connectivity. Correlations are less significant, and sometimes negative, when comparing road network centrality with combined centrality, with the exception of betweenness centrality in the O-D matrix and the spatial network. Mainland cities have less significant correlations than ports on average (around 0.4 compared with 0.9). Closeness and betweenness go in fact in the same positive direction since closeness centrality should be interpreted in the inverse way: higher values mean higher remoteness from the other nodes. This correlation increases from the grid to the matrix and to the spatial network, meaning that road network centrality plays a greater role in intermodalism when the base network takes into account geography, as the maritime network in this dimension is not made of virtual or strait lines but realistic routes avoiding land surfaces. The exclusion of zero values did not alter the results as seen in the figure. Overall, such results are a first proof of the efficiency of our network model, since maritime and land-sea networks are relatively similar with the only difference of including the Australian road network in the former; and given the simplicity and relatively non-hierarchical and circular road network of Australia, it is understandable that road network centrality is not fully correlated with (global) sea-land centrality.

[Insert Table 1 about here] A look at the centrality differentials and their spatial distribution across Australia provides additional clues (Figures 9 and10), with a distinction between mainland cities (Figure 9) and mainland ports (Figure 10). First of all, among all the Australian cities retained in this study, it is generally the case that whatever the network and the indicator considered, the ratio between single and combined centrality is more advantageous for port cities. Inland, non-port cities often exhibit low centrzlity increase or even centrality loss, especially for the southeast region and for betweenness centrality, as most port cities benefit from being directly connected to both networks. However, some port cities in the western region such as Karratha and Geraldton suffer important centrality loss in the land-sea network due to their remoteness and peripherality in the country and especially the road network (see also Alice Springs in the O-D matrix for betweenness). Other examples of disadvantaged port cities are also to be seen in the eastern part, such as Gladstone, Maryborough; Geelong and Launceston in the southern region. Except from the latter city, Tasmanian cities witnessed the largest centrality gains given their islandic situation, cut from the mainland road network and therefore taking full advantage of being connected by water. Surprisingly, no loss occurred for closeness centrality, but important gaps between gaining and losing cities, mostly in the South and among which the largest cities Perth, Melbourne, and Sydney but also peripheral cities such as Darwin and Cairns.

[Insert Figure 9 about here] When it comes to mainland ports as the unit of analysis (many of them being studied in previous Figure 9 as cities), interesting common trends and exception emerge (Figure 10). Not only a higher number of ports could have been analyzed by closeness centrality than by betweenness centrality for technical reasons, we first observe the same "Tasmania effect" with the highest centrality gains for Tasmanian ports to palliate their insularity, this time for both indicators. In terms of betweenness, the gain from adding the maritime network to the road network is mainly observable in the southeast region (except from Perth in the maritime grid), such as for Sydney, Melbourne, and Adelaide, which as large cities still occupy an important situation in the road network. Conversely, poorly connected ports in the road network, such as Darwin, also benefit from this combination, mainly for artificial reasons such as its proximity to Southeast Asian ports and making one of the few direct and shorter links between Australia and its neighbor countries and the hub of Singapore (see [START_REF] Wu | Between the centre and the periphery: The development of port trade in Darwin, Australia[END_REF]. Without commodity flow data on land transport within Australia, our results rely on pure infrastructure that is theoretically and potentially crossable from East to West and North to South. Darwin thus acts as a bridge between the Australian network and the global network, although its gain remains lower than for truly commercial cities in the South.

In terms of closeness centrality, certain large southern ctieis still perform well in most networks, but this time important gains concentrate in the western part of the country, from Port Hedland in the North to Esperance in the South. As an artificial effect of higher grid density along the west coast, such ports were probably advantaged in the maritime grid than in other networks. The low gains of eastern ports on average can also be interpreted as a strength; the inclusion of the maritime network has limited effect given their higher closeness centrality in the road network. In turn, the higher gains of western ports may also reflect their land peripherality so that adding the maritime network greatly improve their intermodal accessibility as a compensation of their remoteness to core urban settlements.

[Insert Figure 10 about here]

Sea-land centrality and the urban hierarchy

We now shift the focus on urban development through the prism of demographic population as a proxy for cities' importance socio-economic importance. We investigate the degree of overlap between centrality measures and urban population (Table 2), distinguishing amongst the host city of port (Pop1) and the nearest larger city (Pop2). One first result is the particular case of island ports (excluding Tasmania), mainly located in the western waters of the country, and which, due to their small size, isolation and absence of road network, exhibit the highest correlation between centrality and population in all networks, but more for Pop1 than for Pop2. Mainland ports have a lower correlation in all networks than island ports, but it remains rather significant (i.e. around 0.4), compared with mainland cities, which have even lower correlations on average. A common result, extremely valuable for responding positively to our initial hypothesis (cf. a combined network is more in line with urban activity than a single network), is that correlations between centrality and population in the sea-land network are always superior, whatever the network considered, to correlations for single networks, should it be maritime (zero correlations) or road (slightly negative correlations). This means that although the modeled road network is not in full accordance with the urban hierarchy measured by population, intermodal centrality is more in line with the urban hierarchy than centrality calculated on either individual networks.

[Insert Table 2 about here]

Conclusion

The application of combined geomatics and graph-theoretical methods to produce a hybrid network made of road infrastructure and maritime flows is a new way forward in both network science and transport studies, notwithstanding noticeable contributions to urban studies. As a first step in that direction, this article focuses on fundamental methodological aspects without which such an effort would not be fully and rigorously justified. The application to Australian port cities provides a number of evidences that are specific to this study area but has the potential to apply to the global sea-land network in future research. However, each territory has its own specificities that ought to be considered and included in the model. In its current state, it has shown interesting results for Australia but still needs some adjustments to fit the specifics of each different road networks found around the globe. Indeed, all countries do not identical road network structures; this can be due to a lack of information but can also be due to socio-economic factors explaining the lack of transportation network's development.

Ongoing applications to other parts of the world already revealed certain issues, such as the heterogeneity of the road hierarchy and categorization that obviate calculations on coastal areas, thus hindering the connection of ports to the road system. For instance, small islands do not necessarily possess motorways or trunk ways but may have several ports connected by lower-rank roads. Another issue is the use of different transport modes in countries such as Vietnam or the Democratic Republic of the Congo. Those two countries have in common to be divided by a large river (i.e. the Mekong and the Congo), so that the network's connectivity is maintained through river crossings via ferries to overcome the lack of bridges along the riverbed. This feature is being implemented in our model with a proper cost weight to reflect the impact of such a diversity and to guarantee the presence of terrestrial links between northern and southern Africa. Lastly, transportation systems are also diverse in terms of their level of network density, from the sparsely populated Russian steppes to the Japanese megalopolis between Tokyo and Fukuoka. The use of the quadtree partitioning method seems to prevent inconsistencies linked to the density of population and transportation infrastructures, but morphologic features such as insularity or, on the contrary, supercontinents such as Afro-Eurasia or the Americas, will inherently bring specific issues. Other improvements related to purely technical aspects, to enhance both optimization and robustness to obtain a fast, deployable model everywhere. For example, coastal areas currently take a significant amount of time to be calculated; this is mainly due to the process of surrounding tiles that, to ensure the seamlessness of those areas, highly increase the operating time and induce a redundancy that must be taken care of.

In the end, the model should be working nicely with OpenStreetMap dataset, which contains all the required features for this kind of network creation [START_REF] Gil | Building a multimodal urban network model using OpenStreetMap data for the analysis of sustainable accessibility[END_REF]. But the idea is to create a version that works as well with other sources such as the Vector Map (VMAP) or any historical data that would enable studies on the period covered by the informations extracted from the Lloyd's List corpus. It will consequently allow us to adapt the model more easily to other road network datasets such as The Vector Map (or VMAP) for the 1980/90 period. As the mechanics and spatial operations will be the same, those adjustments shall be quickly implemented, following an equivalent protocol and modifying only minor features. We also consider weighting the global road network using freely accessible information on the pavement rate at the country level using Central Intelligence Agency World Factbooks 2 . Such an information will allow us to differentiate the score of road network segments depending on the overall quality and navigability of the transport system.

Overall, this article contributes to a somewhat classical issue in transport economics and geographythe sea-land interfacewhile pushing further the technical possibilities of its analysis. It has potential to deepen our understanding of networks in general, such as a bi-modal or multiplex network being both planar (road) and non-planar (shipping), which is still rather uncommon in network science. Computing centrality measures and comparing them with local socio-economic features of the connected nodes will add crucial evidences to the role of network morphology, connectivity, and complexity in the understanding of urban and regional development dynamics. Huge efforts are underway to move towards a more dynamic view of such phenomena, with the extraction of further shipping data and road network data from older sources. As advised by pioneering scholars in the field [START_REF] Rimmer | Maritime Networks. Spatial Structures and Time Dynamics[END_REF]Robinson, 2015), the inclusion of other network types such as railways, airlines, and immaterial linkages (e.g. firms and telecommunications) shall pave the way towards the first-ever analysis of the allencompassing network of networks, to further elucidate the role of specialization and diversification dynamics shaping the evolving world city network. Lastly, another research pathway is the analysis of the maritime centrality of non-port cities situated inland, thanks to their indirect, landside connection with the global shipping network. Whenever possible, road traffic data, available from commodity flow surveys in certain countries, will allow weighting both maritime and landside flows, to complement existing studies focusing only on the continental area [START_REF] Duranton | Roads and trade: Evidence from the U[END_REF].

Table 1 :Figure 1 :

 11 Figure 1: Vessel traffic evolution of Australian ports, 1977-2008 Source: own realization based on Lloyd's List data

Figure 2 :

 2 Figure 2: Vessel traffic distribution among Australian ports, 2008 (Unit: % DWT) Source: own realization based on Lloyd's List data N.B. maximum height = 100%

Figure 3 :Figure 4 :Figure 6 :

 346 Figure 3: Model for the maritime grid construction Source: own realization

Figure 7 :Figure 8 :

 78 Figure 7: Coastal network connection as seen with the case of Brisbane Source: own realization

Figure 9 :

 9 Figure 9: Effects of road network combination on the maritime centrality of Australian cities Source: own realization

Figure 10 :

 10 Figure 10: Effects of road network on the maritime centrality of Australian ports Source: own realization

Containers General cargo Food products Dry bulk Liquid bulk Passengers & Vehicles Other Abbot Point Adelaide Albany(AUS) Ardrossan(AUS) Australia Botany Bay Brisbane Broome Bunbury Bundaberg Burnie Cairns Cape Cuvier Cape Flattery Cossack Pioneer Terminal Dampier Darwin Devonport(AUS) Eden Esperance Exmouth(AUS) Fremantle Geelong Geraldton Gladstone Gove Griffin Terminal Hastings(AUS) Hay Point Hayman Is. Hobart Karumba Kurnell Launceston Legendre Terminal Lucinda Mackay Melbourne Milner Bay Mourilyan Mutineer Terminal Newcastle(AUS) Nganhurra Terminal Port Alma Port Bonython Port Giles Port Hedland Port Hedland Anch. Port Kembla Port Latta Port Lincoln Port Pirie Port Walcott Portland(AUS) Stybarrow Venture Terminal Sydney Thevenard Townsville Useless Loop Wallaroo Weipa Whyalla

Acknowledgements

The research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP/2007-2013) / ERC Grant Agreement n.

[313847] "World Seastems".

Acknowledgements

2 https://www.cia.gov/library/publications/the-world-factbook/

The research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP/2007-2013) / ERC Grant Agreement n.

[313847] "World Seastems".