Maria Choufany

Davide Martinetti

Rachid Senoussi

Cindy E Morris

Samuel Soubeyrand

Spatiotemporal large-scale networks shaped by air mass movements

Keywords: Aerobiology, Air masses dynamics, Connectivity, Spatial network, Spatiotemporal network, Trajectory

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Atmospheric air masses are volumes of air with a defined temperature and water vapor content that have long been known to rule fundamental atmospheric phenomena like weather and air currents. Their composition is mostly inert gases, but both organic and inorganic particles have been found to linger in high-altitude air as a consequence of the constant interaction of air masses with the earth's surface below Email addresses: maria.choufany@inra.fr (Maria Choufany), davide.martinetti@inra.fr (Davide Martinetti), rachid.senoussi@inra.fr (Rachid Senoussi), cindy.morris@inra.fr (Cindy E. Morris), samuel.soubeyrand@inra.fr (Samuel Soubeyrand) them. A non-exhaustive list includes gases and minerals like wildfire smoke, radioactive material, dust, sand, volcanic ash and sea salt, but also living organisms such as pollen, fungal spores, bacteria, virus and small insects. Despite the relative sparse density of these particles with respect to the volume of an air mass, their presence and transportation across the planet has proven to have strong effects on many phenomena impacting human health and safety (pollen [START_REF] Mahura | Elevated birch pollen episodes in Denmark: contributions from remote sources[END_REF][START_REF] Šauliene | Application of backward air mass trajectory analysis in evaluating airborne pollen dispersion[END_REF][START_REF] Bogawski | Detecting distant sources of airborne pollen for Poland: Integrating back-trajectory and dispersion modelling with a satellite-based phenology[END_REF], dust concentrations [START_REF] Khaniabadi | Impact of Middle Eastern Dust storms on human health[END_REF][START_REF] Aciego | Dust outpaces bedrock in nutrient supply to montane forest ecosystems[END_REF], nuclear byproducts [START_REF] Moroz | Predictions of dispersion and deposition of fallout from nuclear testing using the NOAA-HYSPLIT meteorological model[END_REF][START_REF] Rolph | Modeling the fallout from stabilized nuclear clouds using the HYSPLIT atmospheric dispersion model[END_REF], human, animal and plant epidemics [START_REF] Leyronas | Assessing the aerial interconnectivity of distant reservoirs of sclerotinia sclerotiorum[END_REF][START_REF] Wang | Long-distance spore transport of wheat stripe rust pathogen from Sichuan, Yunnan, and Guizhou in southwestern China[END_REF][START_REF] Aylor | The role of intermittent wind in the dispersal of fungal pathogens[END_REF][START_REF] Mundt | Aerial dispersal and multiple-scale spread of epidemic disease[END_REF][START_REF] Sadyś | Back-trajectories show export of airborne fungal spores (Ganoderma sp.) from forests to agricultural and urban areas in England[END_REF][START_REF] Hiraoka | Seasonal analysis of microbial communities in precipitation in the greater Tokyo area Japan[END_REF], air pollution (Liu et al., 2018b,a;[START_REF] Talbi | Assessment of annual air pollution levels with PM1, PM2.5, PM10 and associated heavy metals in Algiers, Algeria[END_REF], and rainfall [START_REF] Chen | Analysis of paths and sources of moisture for the south China rainfall during the presummer rainy season of 1979-2014[END_REF][START_REF] Armon | Synoptic-scale control over modern rainfall and flood patterns in the levant drylands with implications for past climates[END_REF][START_REF] Rabinowitz | An investigation of atmospheric rivers impacting heavy rainfall events in the North-Central Mississippi River Valley[END_REF]).

The rise in the number of publications on these subjects suggests a growing interest of the scientific community on the effects of air-mass movements on the biosphere, that has surely been boosted by recent available developments, such as the Hybrid Single-Particle Lagrangian Integrated Trajectory model (HYSPLIT, [START_REF] Stein | NOAA's HYSPLIT atmospheric transport and dispersion modeling system[END_REF]), allowing reconstruction of actual air-mass movements at rather fine geographical and temporal scales and with a global cover.

The vast majority of studies focused on isolated events, such as dust storms or peaks of air pollutants, that are rather concentrated in time (from few hours to few weeks) and/or space (just a few locations such as cities). Nonetheless, the movement of air masses is expected to have impacts on a broader spatiotemporal scale, as reviewed in recent studies [START_REF] Leyronas | Assessing the aerial interconnectivity of distant reservoirs of sclerotinia sclerotiorum[END_REF][START_REF] Margosian | Connectivity of the American agricultural landscape: assessing the national risk of crop pest and disease spread[END_REF]. The purpose of the present paper is then to propose a mathematical framework for studying air-mass movements on large spatiotemporal scales, under the hypothesis that these movements can create stable and recurrent connections between distant portions of a territory. The very nature of these connections will be further specified throughout the manuscript, but as a general rule we will consider that any pair of points (or areas) in space can have a certain degree of connection, regardless of their geographic distance, provided that there are recurrent air-mass trajectories that connect the two points (or areas). The direction and strength of these connections will be estimated by looking at the trajectories linking every pair of points/areas and weighting them according to appropriate measures. In this perspective, it seems natural to resort to graph and network theory, since the formalism of nodes and edges provides an adequate environment for describing complex connections and can further be used to deepen into the topology of the constructed networks in order to infer interesting properties of the graphs, such as the presence of hubs.

From a generic statistical point of view, we aim to (i) estimate the weighted and directed edges of a graph using a sample of trajectories of individuals traveling through the space formed by the nodes of the graph, and (ii) characterize the estimated graph based on relevant statistics.

In the following sections, we first introduce the definitions and properties that will allow us to describe and then estimate connections between points/areas in space via spatiotemporal trajectories. Then, we propose several types of measures to model diverse types of connections. The expected output consists of a spatiotemporal graph describing the network of links induced by trajectories. It's worth noting that our approach is meant to infer connectivity induced by air-mass movements and it is readily applicable to HYSPLIT-type data, but we have maintained a sufficient level of generality to be applied to other phenomena, provided that trajectory data are available (e.g. animal trajectories). Finally, we apply our method to two case studies concerning the coastline of the Mediterranean sea and the French region of Provence-Alpes-Côte d'Azur. The two case studies have different spatiotemporal granularities and they will be used to provide examples of application of the proposed methodology.

Framework for the definition of trajectory-based networks

In this section we show how a set of trajectories evolving within space during a finite time interval can be used to construct pertinent spatiotemporal networks. We first recall some basic definitions related to networks (Section 2.1) and then propose a statistical methodology to infer the network structure from a data set of trajectories (Section 2).

Network theory

Network theory (a.k.a. graph theory) is a mathematical formalism introduced by Leonhard Euler to describe the famous Königsberg bridge problem [START_REF] Newman | The structure and function of complex networks[END_REF][START_REF] Strogatz | Exploring complex networks[END_REF][START_REF] West | Introduction to graph theory[END_REF]. The two basic components of a network are a set of nodes linked by a set of edges. Nodes can represent a variety of things, such as persons, regions, computers, neurons, etc., while edges are used to describe the connections between those nodes. Formally, a network G " pV, Eq is defined as a set of nodes (or vertices)

V " tv 1 , v 2 , . . . , v N u connected by a set of edges E " te ij u i,jPt1,...,N u . A natural way of representing a network is given by means of a N ˆN square matrix M , usually referred to as an adjacency matrix, whose term pi, jq, M ij , is non-zero an edge exists between i and j. By convention, adjacency matrices are defined to have an empty diagonal (i.e. M ii " 0, i P t1, . . . , N u), meaning that nodes cannot be self connected. If M is symmetrical (i.e. M ij " M ji , i, j P t1, . . . , N u), then the network is said to be undirected, and directed otherwise. If M ij P t0, 1u, the network is said to be binary, meaning that an edge between two nodes i and j either exists or does not. Otherwise, if M ij P R, the network is said to be weighted, meaning that the edge between nodes i and j are more or less connected.

In this paper, a network is said to be spatial [START_REF] Barthélemy | Spatial Networks[END_REF] when nodes correspond to geographic locations, while we use the term temporal [START_REF] Holme | Temporal networks[END_REF] to refer to networks where edge values can change over time. Finally, we will use the term spatiotemporal network to refer to network that are simultaneously spatial and temporal, under the constraint that nodes cannot change position, neither appear nor disappear over time. The networks considered in this paper also fall into the rather generic definition of spatiotemporal networks. If the spatial qualifier means that the nodes of the networks represent fixed geographical locations, the temporal qualifier is more complex. Indeed, temporal networks are generally divided in the literature into two main classes, namely contact graphs or interval graphs [START_REF] Holme | Temporal networks[END_REF]. The former type refers to networks where edges represents instantaneous contacts between nodes (Figure 1 and B in panel (b) is active during the time intervals (0,3) and (6,8). For contact networks (c), the edges are quantitatively more or less active across time, and the quantity of activity of any edge is described by a temporal function.

Flows and trajectory segments

We consider a function Φ : R ˆR ˆΩ Ñ Ω, usually called flow on the spatial domain Ω of R d , satisfying the following properties:

$ & % Φpt, s, xq " Φpt, t 1 , Φpt 1 , s, xqq Φ ´1pt, s, .q " Φps, t, .q, ( 1
)
where s, t, t 1 P R and x P Ω. For fixed t and s, the flow Φpt, s, ¨q is a spatial transformation. For fixed x and varying s or t, the function gives a forward or backward trajectory of a particle over Ω between times t ^s " infpt, sq and t _ s " suppt, sq. If s ď t, y " Φpt, s, xq gives the future location at time t of the particle presently located at x at time s. Contrarily, if s ě t, y " Φpt, s, xq gives the location at past time t that was occupied by the particle located at x at present time s. Φpt, s, .q is assumed to be a bijective mapping meaning that particles following distinct trajectories cannot be at the same location at the same time.

In general, a flow is defined with respect to a possibly time-dependent vector field F over R ˆΩ, as the solution u : R Ñ Ω of an ordinary differential equation (see e.g. Hamilton's equations in classical mechanics) with specified initial condition at a specified time s:

$ & % du dt ptq " F pt, uptqq upsq " x, ( 2
)
where F is continuous and Lipschitzian over R ˆΩ. In the setting introduced above, Φpt, s, xq " uptq with Φps, s, xq=u(s)=x. The solution u represents the trajectory of the particle located at x at time s.

Varying the initial condition in System (2), i.e. varying s and x, leads to consider pieces of trajectories of all particles which dynamics are governed by the vector field F . In this article, the vector field F will not be made explicit, but we will consider samples of trajectory segments (defined below) for constructing trajectory-based networks.

Definition 2.1. The trajectory segment associated to the flow Φ over the time interval ∆ ts " rt ^s, t _ ss, s, t P R, for a particle located at x P Ω at time s is defined as follows:

Γpt, s, xq " tpt 1 , Φpt 1 , s, xqq : t 1 P ∆ ts u.

(3)

If s ă t (resp. s ą t), Γpt, s, xq is a forward (resp. backward) trajectory segment. In this article, we are mainly interested in backward trajectories, but the framework presented here encompasses forward trajectories as well.

Example 2.1. The notions of flow and trajectory segment can be adapted to cope with air mass trajectories over the Earth surface. In this case, the spatial domain Ω representing the Earth surface is the sphere S 2 in R 3 . If in addition, air masses are characterized by altitude and temperature evolving in space and time, then Ω " S 2 ˆR`ˆR , where R `(resp. R) is the domain of the altitude (resp. temperature) coordinate.

Example 2.2. Animal movements and behaviour activities can also be represented with the notions of flows and trajectory segments, providing, for instance, the animal locations and the covariate value indicating whether animals are feeding or not. In this case, Ω " R 2 ˆt0, 1u, where 1 stands for 'the animal is feeding' and 0 otherwise. The use of a binary variable for describing the feeding activity may require the use of stochastic processes or generalized functions undergoing dynamic analog to the System

(2) for constructing the flow if it is defined with respect to a vector field F .

Pointwise and integrated connectivities

Trajectory-based networks are grounded on the notion of connectivity used as a quantitative, directed measurement of edges between graph nodes. In this aim, we first define the pointwise connectivity as a measure (or submeasure), in the mathematical sense, of the connectivity between a subset A and a point

x of Ω induced by the trajectory segments Γpt, s, xq of a particle located at x at time s. Then, we use the pointwise connectivity to define the integrated connectivity between two subsets A and B of Ω over a temporal domain ∆ of R (∆ can be the union of disjoint intervals).

A Þ Ñ Ψ C pA | t, s, xq is only a submeasure on Ω since Ψ C pA Y A 1 | t, s, xq ď Ψ C pA | t, s, xq `ΨC pA 1 | t, s, xq
for disjoint sets A and A 1 of BpΩq.

Remark 1. This example based on the simple contact between sets can be considered as too strict from a statistical and measure-theory perspective since the length or the duration of a contact may be null.

Instead, a positive constraint on contact length for example can be used to define another version of the contact-based pointwise connectivity: Equation (4) could then be replaced by

Ψ C pA | t, s, xq " 1 tLpAtsXΓpt,s,xqqą0u ,
where LpA ts X Γpt, s, xqq denotes the length of the curve Γ within A. The length operator L will be made explicit in Example 2.5.

Example 2.4. The duration-based pointwise connectivity is defined by:

Ψ D pA | t, s, xq " ż ∆ts 1 tΦpv,s,xqPAu dv, ( 5
)
to measure the duration spent by the particle in A during ∆ ts .

Example 2.5. The length-based pointwise connectivity is defined by:

Ψ L pA | t, s, xq " ż ∆ts 1 tΦpv,s,xqPAu ||∇ v Φpv, s, xq||dv, (6)
where ∇ v Φpv, s, xq stands for the gradient of the flow Φ with respect to the time variable v and || ¨|| denotes the Euclidean norm. Ψ L pA | t, s, xq measures the distance travelled within A by the particle during ∆ ts .

Example 2.6. The pointwise connectivity based on local volume is defined by:

Ψ V pA | t, s, xq " ż ∆ts 1 pΦpv,s,xqPAq | detpJ x Φ pv, sq|dv (7)
where detpJ x Φ pv, sqq is the determinant of the Jacobian matrix (with respect to x) of the spatial transformation Φpv, s, ¨q. The absolute value | detpJ x Φ pv, sq| of the Jacobian determinant at x gives the ratio by which the function Φpv, s, ¨q expands/shrinks infinitesimal volumes around location x into infinitesimal volumes around location Φpv, s, xq.

In other words, Ψ V pA | t, s, xq assesses how particle density increases or decreases from x to A along the time interval ∆ ts . Intuitively, if n particles are initially in A and if the infinitesimal volume around any of these particles tends to shrink from A to x, then one expects a high concentration of particles in a fixed volume around x and, therefore, a high connectivity from A to x. Conversely, if the infinitesimal volume around a particle tends to expand from A to x, then one expects a lower concentration of particles in the same fixed volume around x and, therefore, a lower connectivity from A to x.

More sophisticated specifications of the pointwise connectivity can be proposed by incorporating spatio-temporal covariates in its formulation, like in the following examples.

Example 2.7. Let G denote a time-varying vector field defined over R ˆΩ. The pointwise connectivity based on the external vector field G is defined by:

Ψ G pA | t, s, xq " ż ∆ts 1 pΦpv,s,xqPAq | ă ∇ v Φpv, s, xq, Gpv, Φpv, s, xqq ą | dv
where ă ∇ v Φpv, s, xq, Gpv, Φpv, s, xqq ą is the scalar product between the gradient with respect to the time variable v of the flow Φ and the vector field G. Larger the average collinearity in A between the instantaneous movement of the particle and the simultaneous direction of the vector field G, higher the connectivity between A and x. For instance, if Φ gives the movement of air masses and G provides the intensity and the direction of a continuous release of specific particles, then the connectivity will be high (resp. low) if the movement of the air in A and the movement of particles released in A are approximately collinear (resp. orthogonal).

Example 2.8. Let Z and Z be positive real valued spatio-temporal functions defined over R ˆΩ. The pointwise connectivity based on Z and Z is defined by:

Ψ Z, Z pA | t, s, xq " Zps, xq ż ∆ts 1 pΦpv,s,xqPAq Zpv, Φpv, s, xqqdv. ( 8
)
This form of pointwise connectivity may represent, for example, (i) the negative effect of the altitude of the air mass when it is above A on the recruitment of specific particles from the ground, and (ii) the positive effect of rainfall at ps, xq on the deposition of particles from the air mass to the ground (see Figure 3). Thus, lower the average altitude of the air mass above A and more intense the rainfall at ps, xq, larger the contribution to the connectivity between A and x. This is expressed in Equation (8)

as follows: (i) Z is defined as the binary function indicating whether or not the altitude of the air mass 6) to account for a supplementary effect of the distance travelled within A on the pointwise connectivity.

(
Each pointwise connectivity defined above can be used for defining the integrated connectivity, which measures the quantitative directional link between two subsets A and B of BpΩq generated by trajectories of particles located in B at times belonging to the temporal domain T .

Definition 2.3. Let A and B be two sets of BpΩq and T a subset of the temporal domain R. The δ-lag integrated connectivity linking B to A over T is defined by:

Ψ p2q ν,δ pB ˆA | T q " ż T ˆB ΨpA | s `δ, s, xqνpds, dxq, (10)
where δ P R and ν is a measure on R ˆΩ. ν is both discrete in space and time with a mass only at a countable collection of space-time points.

Trajectory-based network

Definition 2.4. A trajectory-based network generated by Ψ p2q ν,δ (given by Equation ( 10)) over the temporal domain T Ă R, is a graph whose nodes A i , i " t1, .., Iu, are disjoint sets of Ω in BpΩq and whose directed edges are weighted by integrated connectivities M ij " Ψ p2q ν,δ pA i ˆAj | T q, 1 ď i, j ď I and i ‰ j.

Definition 2.4 corresponds to a spatial trajectory-based network evaluated over the fixed temporal domain T . It can be extended in different ways to obtain spatiotemporal analogs. For example, if T 1 , . . . , T K denote K disjoint but successive time intervals with equal lengths, then the sequence of trajectory-based networks generated by Ψ p2q ν,δ p¨ˆ¨| T k q, k " 1, . . . , K, forms a spatiotemporal trajectorybased network that can be analyzed to assess how connectivities across space are changing with time. This is one of the issues considered in Section 4.2. 1.0 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q (a)

x [1 ] x [2 ]
Ti m e q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q x[1] 1.0 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q (a)

q q q q q q x[1] x[2] q q q q q q x[1] x[2] Time (b)

Estimation of integrated connectivities

In practice, the integral defining the integrated connectivities between subsets of Ω (Definition 2.3) cannot be analytically computed in general, but can be estimated from a sample of trajectories. For instance, the integrated connectivity Ψ p2q ν,δ pB ˆA | T q can be estimated by its empirical counterpart obtained by importance sampling, say p

Ψ p2q ν,δ pB ˆA | T q: p Ψ p2q ν,δ pB ˆA | T q " 1 N 1 N 1 N ÿ k"1 N 1 ÿ l"1 ΨpA | s k `δ, s k , b l q, ( 11
)
where s 1 , . . . , s N P T and b 1 , . . . , b N 1 P B are times and locations, respectively, randomly drawn under the measure ν restricted to T ˆB, and ΨpA | s k `δ, s k , b l q is the pointwise connectivity associated to the trajectory of the particle located at b l at time s k and observed over ∆

s k `δ,s k " rs k ^sk `δ, s k _ s k `δs.
If ν is constant, other classical numerical approaches can be applied to approximate the integral, such as an hybrid approach in which the mid-point rule is applied in time and a regular point process is used in space. In such a case, the integrated connectivity estimator is also given by Equation (11).

Example 3.1. Using Equation ( 11), the contact-based pointwise connectivity in Example 2.3 is estimated by:

p Ψ p2q C,δ pB ˆA | T q " 1 N 1 N 1 N ÿ k"1 N 1 ÿ l"1 1 tA s k `δ,s k XΓps k `δ,s k ,b l q‰Hu , ( 12
)
where A s k `δ,s k " ∆ s k `δ,s k ˆA. Thus, p Ψ p2q C,δ pB ˆA | T q is simply the proportion of sampled trajectories intersecting A.

Example 3.2. Using Equation ( 11), the duration-based pointwise connectivity in Example 2.4 is estimated simply by the average duration of the intersections between the sampled trajectories and A.

Example 3.3. Using Equation ( 11), the length-based pointwise connectivity in Example 2.5 is estimated simply by the average length of the intersections between the sampled trajectories and A.

Applications

In this section, we applied our general framework to the flow of air mass movements. Indeed, these movements compiled over years were used to characterize climatic patterns [START_REF] Hondula | A back-trajectory and air mass climatology for the Northern Shenandoah Valley, USA[END_REF] and to describe the transport of pollutants [START_REF] Pérez | Applications of air mass trajectories[END_REF]. We show now how to deploy our approach for constructing air-mass movement networks in two real geographical contexts, namely the coastline of the Mediterranean Sea and the French region of Provence-Alpes-Côte d'Azur. These two examples have been chosen in order to prove the flexibility of our approach to different situations and geographical scales.

Case study regions and network construction

The first study region corresponds to the coast of the Mediterranean Sea, ranging approximately 1,600 km from north to south and 3,860 km from east to west. The temperate climate of the chosen region is strongly influenced by the presence of the Mediterranean Sea, with mild winters, hot summers and relatively scarce precipitations events. The landscape is characterized by coastal vegetation, typically shrubs and pines, and densely populated areas with intensive crop production of wheat, barley, vegetables and fruits, especially olive, grapes and citrus. In this paper, we characterize recurrent movements of air masses through the Mediterranean region by defining a grid with mesh size 74 km covering the coastline from 5 km up to 250 km inland from the coast, including the four largest islands (namely Sicily, Sardinia, Cyprus and Corsica). Thus, we divide the region into 604 cells, where the centroids of the cells will be used as arrival locations of air-mass trajectories and will correspond to the nodes of the constructed network.

The second study region corresponds to the French region of Provence-Alpes-Côte d'Azur (PACA, hereafter), located in the south-eastern part of France and characterized by a rather complex landscape formed by a densely-populated coastline, agricultural lands (high-value-crops with fruit and olive orchards, vineyards, vegetable cultivation and horticulture), and natural mostly-alpines areas. The choice of this particular region is justified in the context of a research project aimed at assessing the potential longdistance dissemination of phytopathogenic bacterial populations that are known to be transported by air currents. The bacteria of interest (e.g., Pseudomonas syringae) can be lifted in to the air from a source location and then be passively transported by air masses until they are deposited back to land onto a different, far away sink location. Since the life cycles of the considered species of bacteria are strongly linked to the water cycle [START_REF] Morris | The life history of the plant pathogen Pseudomonas syringae is linked to the water cycle[END_REF], we naturally partitioned the study area in a way that fit this assumption. Hence, we considered the 294 watersheds of the PACA region to define the sites that will later constitute the nodes of the constructed network. Since watersheds have irregular shapes and varying sizes, we selected a certain number of arrival locations per watershed (between 1 and 10 and proportionally to the watershed area) in order to cover the watersheds consistently according to the relative importance of their size and estimate the integrated connectivities. In total, a set of 833 arrival locations for air-mass trajectories was generated.

Once the arrival points for the two study regions have been established, we turned to the computation of air-mass trajectories arriving at the prescribed locations using the Hybrid Single-Particle Lagrangian Integrated Trajectory model (HYSPLIT, [START_REF] Stein | NOAA's HYSPLIT atmospheric transport and dispersion modeling system[END_REF]). The HYSPLIT model can be fed with meteorological data from the Global Data Assimilation System files with a 0.5-degree spatial resolution (GDAS 1 ) and was tuned by us to return 48-hours backward air-mass trajectories arriving at the prescribed locations at an altitude of 500 m above mean sea level. A single trajectory consists of a vector containing the hourly positions (longitude, latitude and altitude) visited by the air mass before arriving at the specified location and time. Air-mass trajectories have been computed for every arrival location (604 for the Mediterranean region and 833 for the PACA region) and for every day between January 1, 2011 and December 31, 2017 (arrival hour is 12:00 GMT). The total number of computed trajectories is 1,543,220 for the Mediterranean region and 2,128,315 for the PACA region.

The final step for the construction of the networks is the estimation of the adjacency matrices of the networks, based on the methodology presented in the previous sections. To do that, for each pair of subsets of the spatial domain, we used the daily 48-hours backward trajectories arriving at the locations sampled within the receptor subset, and computed the contact-based estimator (see Example 3.1). The subsets of the spatial domain are the watersheds for PACA and circular buffers of radius 20 km for the Mediterranean region, as in [START_REF] Leyronas | Assessing the aerial interconnectivity of distant reservoirs of sclerotinia sclerotiorum[END_REF] In this work we will consider networks corresponding to three temporal contexts: (i) the spatial networks obtained when T is the entire period 2011-2017, (ii) the yearly spatiotemporal networks formed by the seven spatial networks obtained when T 1 encompasses the year 2011, T 2 encompasses 2012 and so on, and (iii) monthly spatiotemporal networks formed by the twelve spatial networks obtained when T 1 represents every January from 2011 to 2017, T 2 every February from 2011 to 2017, and so on.

Network analysis

The constructed networks are directed and weighted by contact-based connectivities generated by air mass trajectories. They are inherently complex by the sheer amount of spatial and temporal information that they encompass. Hence, there is no easy way of representing the results either graphically or numerically, without compromising the original complexity of the networks. While a comprehensive physical study of the spatiotemporal properties of these networks goes beyond the scope of the paper, we explore the estimated trajectory-based networks by looking at some generic properties through the 1 https://www.ncdc.noaa.gov/data-access/model-data/model-datasets/global-data-assimilation-system-gdas following indices:

• Diameter: the longest of all possible shortest paths between any two pair of nodes.

• Density: the ratio between the sum of all edge weights and the number of all possible edges [START_REF] Liu | Complex discovery from weighted ppi networks[END_REF].

• Transitivity (also known as clustering): the equivalent definition of density, but applied to triplets of nodes instead of pairs of nodes [START_REF] Opsahl | Clustering in weighted networks[END_REF].

• Shortest path: characterized by the average and standard deviation of the computed shortest path between any possible pair of different nodes [START_REF] Newman | Scientific collaboration networks. ii. shortest paths, weighted networks, and centrality[END_REF].

• Small worldness: refers to the property of a network of being highly clustered and having relatively short shortest paths. It is computed as the ratio between the normalized clustering and the normalized average shortest path distance [START_REF] Li | The structure of weighted small-world networks[END_REF][START_REF] Colon-Perez | Small worldness in dense and weighted connectomes[END_REF].

• Scale-free property: The degree of a node in terms of the total number of edges entering and exiting from it, and for directed networks it can be decomposed in the incoming and outgoing degree, respectively. The degree distribution is the empirical distribution of the degree of a network and it said to be scale free when it approximately follows a power law distribution, i.e. P pkq " k p´αq , where P pkq represents the probability of a node having degree equal to k [START_REF] Barabási | Scale-free networks[END_REF][START_REF] Barabási | Emergence of scaling in random networks[END_REF]. Some authors impose that the α parameter of the power law distribution has to fall within the interval r2, 3s [START_REF] Barabási | Network science[END_REF]. Thus, a network is scale free when most of its nodes have low degree, while the probability of having nodes with very high degree is not negligible (fat right tail of the distribution). Nodes with very high degree play a crucial role in dynamics conditional on networks and are often referred as hubs [START_REF] Liu | Controllability of complex networks[END_REF].

• Degree correlation: in directed networks, it accounts for the correlation between the incoming and the outgoing degree of a node. Networks with positive (resp. negative) degree correlation foster (resp. hamper) epidemic spread [START_REF] Pautasso | Disease spread in small-size directed trade networks: the role of hierarchical categories[END_REF]. The two spatial trajectory-based networks representing the strength of tropospheric connections in the Mediterranean region and PACA during the entire period 2011 to 2017 are represented in Figure 4. In order to highlight those edges that represent strong connections, we depicted them with darker shades of color, while we did not draw the connections that had a weight of less than 0.3. It can be seen that the strongest connections tend to link nodes that are geographically close, but nonetheless moderate connections also exist between rather distant nodes. This is confirmed by small values of the average shortest path distance 8.20 ˆ10 ´4 for the Mediterranean and 2.57 ˆ10 ´4 for PACA, and high values of the transitivity index (0.74 for the Mediterranean region and 0.99 for PACA), as shown in the first lines of Tables 1 and2. An interesting difference between the two networks is that the one for PACA has a very negative degree correlation (´0.85), meaning that nodes having a high incoming degree will have low outgoing degree, and vice versa. On the other hand, for the Mediterranean network, the value of the degree correlation is moderately positive (0.31), meaning that nodes having high a incoming degree tend to have also high a outgoing degree. Qualitatively, the indices provided in Tables 1 and2 are overall more variable for the monthly spatiotemporal trajectory-based networks than for the yearly ones. Thus, focusing in what follows the monthly network, we investigate possible seasonal patterns by using the complete-linkage hierarchical clustering method [START_REF] Ferreira | A comparison of hierarchical methods for clustering functional data[END_REF]. We applied the clustering using the Euclidean distance over the 8-dimensional space formed by the 8 indices provided in Tables 1 and2.

Results

Mediterranean region

Diam

For the Mediterranean region, the dendogram in Figure 5 For PACA, the dendogram in Figure 7(a) can be used to identify three distinctive periods: summer months (from June to August), winter and spring months (from December to April, plus September that can be considered as an outlier from a chronological viewpoint) and a set of transition months between the two previous periods (May, October and November). incorporate the uncertainty about the trajectory reconstruction impacted by an eventual preliminary step consisting in estimating the parameters of the above-mentioned probabilistic model.

Concerning the application treated in this article, we observed distinct seasonal patterns in the temporal variation of the networks covering the Mediterranean coastline and PACA. In the former case, the networks corresponding to the four clusters shown in Figure 5(b-e) exhibit different spatial patterns of hubs (in terms of location and size) and different trends in the main connectivity directions. In the latter case, the differences between the three networks identified with the clustering approach are mostly related to connectivity amplitude. It would be interesting to explore whether this observation made at two very different spatial scales and resolutions generally holds by studying regions of size similar to PACA all along the Mediterranean coastline.

In the long-term context of our applied research projects connected to aerobiology, the construction and exploration of networks generated by air-mass movements are a way to unravel epidemiological dynamics (and the resulting genetic patterns) of microbial pathogens disseminated at long distance via air movements in the troposphere (see [START_REF] Leyronas | Assessing the aerial interconnectivity of distant reservoirs of sclerotinia sclerotiorum[END_REF], for a proof of concept). Indeed, even if the pathogen is not explicitly taken into account by the framework proposed in this article, the description of connectivities that it offers provides us a proxy of airborne pathogen movements over long temporal terms and large spatial scales. This proxy is a mean to understand pathogen transportation and to anticipate its long distance dissemination. Specifically, network indices such as those calculated in this article can be associated with particular epidemiological properties such as the probability of long-distance transport of pathogens [START_REF] Moslonka-Lefebvre | Networks in plant epidemiology: from genes to landscapes, countries, and continents[END_REF][START_REF] Jeger | Modelling disease spread and control in networks: implications for plant sciences[END_REF][START_REF] Pautasso | Network epidemiology and plant trade networks[END_REF]. For instance, for plant pathogens, recent studies [START_REF] Nicolaisen | Fungal communities including plant pathogens in near surface air are similar across Northwestern Europe[END_REF][START_REF] Bowers | Seasonal variability in bacterial and fungal diversity of the near-surface atmosphere[END_REF][START_REF] Aho | Spatiotemporal patterns of microbial composition and diversity in precipitation[END_REF] showed that airborne populations of bacteria and fungi are rather constant across the years, while higher diversity can be observed in different seasons. This statement resonates with our analyses where we observed clear seasonal signals in the estimated monthly spatiotemporal networks in Section 4.3 whereas the yearly signals were less obvious.

Finally, the networks estimated using our approach could be a basis for developing epidemiological models (explicitly handling the pathogen) incorporating long-distance dissemination conditional on recurrent air-mass movements. Such models could be exploited to set up surveillance strategies for early warning and epidemic anticipation in order to help reduce the impacts of airborne pathogens on human health, agricultural production and ecosystem functioning [START_REF] Mundt | Aerial dispersal and multiple-scale spread of epidemic disease[END_REF].

Figure 1 :

 1 Figure1: Types of temporal networks. The time of activation is indicated within the grey bar next to the edges (ranging between 0 and 8). For contact networks (a), edges activate only for one instant at the time and are marked with black vertical lines inside the grey bars. For example, in panel (a), the edge between nodes A and B is only active at instants 2, 4, 6 and 7. For interval networks (b), edges can be activated during an interval of time. For example, the edge between A

  Definition 2.3 encompasses connectivities generated by either forward or backward trajectories, depending on the sign of δ. The use of a unique duration |δ| could be relaxed to account for space-time heterogeneities in the duration of trajectories. It could even be infinite by introducing a measure over time like in Equation (9).The measure ν in Definition 2.3 can be continuous, discrete or hybrid over R ˆΩ. Indeed, if particles of interest are air masses, then B can be considered as continuously filled in space and time. Conversely, if particles of interest are animals of a specific species, then animals occupy only punctual locations in B at each time and the measure x Þ Ñ νpds, dxq, given s, is discrete in Ω, whereas the temporal component of ν is continuous. Another examples occurs when the time s corresponds to death times of animals, then

Figure 2 :

 2 Figure 2: Illustration of contact-based, duration-based and length-based pointwise connectivities (resp. Ψ C , Ψ D and Ψ L ) between the elliptic spatial domain A Ă R 2 and different spatial points x at time s " 1, for ∆ts " r0, 1s. The left curve on panel (a) never enters the domain A. The middle curve on panel (a) enters A (red part of the curve) over a relatively long duration (as shown by panel (a)) but a short distance (as shown by panel (b)). The right curve on panel (a) enters A over a shorter duration but a longer distance. Thus, Ψ C pA | t, s, xq, Ψ D pA | t, s, xq and Ψ L pA | t, s, xq are zero for the left curve; Ψ C pA | t, s, xq " 1 for the two other curves; Ψ D pA | t, s, xq is larger for the middle curve than for the right one, whereas Ψ L pA | t, s, xq is larger for the right curve than for the middle one.

Figure 3 :

 3 Figure 3: Illustration of pointwise connectivity based on a covariate measured along the trajectory (see Example 2.8). In this illustration, the passage of the particle in the elliptic spatial domain A contributes to the pointwise connectivity (red part of the curve in panel (a)) only when the particle is at an altitude lower than a threshold value (grey part of the curve in panel (b)).

Figure 4 :

 4 Figure 4: Networks weighted by contact-based connectivities generated by air mass trajectories between (a) the 604 sampled circular areas around the Mediterranean sea and (b) the 294 watersheds of PACA. Edges with weights lower than 0.3 are not drawn.

  Figure 6(a), summer months are characterized by high diameter and density, low values of transitivity, small-worldness and the lowest values of degree correlation. Winter months have lower diameter and show high values of small-worldness due to its high values of clustering and low values of average shortest path distances. Fall months have lower values of density and small-worldness due to its low values of clustering and high values of average shortest path distances. Finally, the group of transition months show high values of density and degree correlation.

  Figure 7(b-d) illustrates the differences between the networks derived from this clustering. The summer network is largely more connected than the two other networks, and the transition months, surprisingly, do not lead to intermediate connectivities but to the lowest connectivities. Based on Figure 6(b), the group of summer months is characterized by a high diameter, density and average shortest path, low values of transitivity, small-worldness and the lowest degree correlations (yet still negative). Winter and spring months have significantly lower diameter, density and average shortest path distances. Finally, the group of transition months show the highest values of small-worldness, due to their high values of clustering and low values of average shortest path distances.

Figure 5 :

 5 Figure 5: (a): Dendrogram of the months obtained from a hierarchical cluster analysis of the Mediterranean spatiotemporal network based on the monthly dissimilarities of the indices presented in Table 1. (b), (c), (d) and (e): Networks corresponding to the four identified clusters where one displayed only the edges between the nodes connected more than 10 days per month via the air mass trajectories.

Figure 6 :

 6 Figure 6: Boxplot for the different indices (Diameter, density, transitivity, shortest path (mean), small worldness, degree correlation) obtained from (a) the four clusters identified for the Mediterranean region (see Figure 5) and (b) the three clusters for PACA (see Figure 7).

Figure 7 :

 7 Figure 7: (a): Dendrogram of the months obtained from a hierarchical cluster analysis of the PACA spatio-temporal network based on the monthly dissimilarities of the indices presented in Table 2. (b), (c) and (d): Networks corresponding to the three identified clusters where one displayed only the edges between the nodes connected more than 10 days per month via the air mass trajectories.

  Diverse types of the pointwise connectivity can be constructed, either using trajectory segments generated by Φ, or directly using Φ. Specific pointwise connectivities can include environmental covariates and even covariates associated to very the movements of particles. Below, we give several examples of such specifications. Some of these examples are graphically represented in Figure2. Most examples are particularly relevant when Ω is a simple geographic domain and when Φ defines movements of individuals (e.g., air masses, animals or particles) within Ω.

Definition 2.2. Let x P Ω and A P BpΩq, where BpΩq is a σ-algebra of subsets of Ω. The pointwise connectivity associated to the flow Φ is defined as a real valued function Ψ on BpΩq ˆR ˆR ˆΩ, conveniently denoted by ΨpA | t, s, xq, where A Þ Ñ ΨpA | t, s, xq is a measure or a submeasure on Ω for each t, s, x Example 2.3. The contact-based pointwise connectivity is defined by: Ψ C pA | t, s, xq " 1 tAtsXΓpt,s,xq‰Hu , (4) where A ts " ∆ ts ˆA and 1 denotes the indicator function. Ψ C pA | t, s, xq indicates whether or not the particle whose movement in Ω is governed by Φp¨, s, xq hit A during the time interval ∆ ts . Note that

  If in Example 2.8, the altitude of the air mass is incorporated as the third coordinate of Φ and A is a 3D-domain vertically limited by the threshold value h, then, Equation (8) is simplify reduced to Equation (5). In the same vein, Example 2.8 can also be modified by adding within the integral the term ||∇ v Φpv, s, xq|| arising in Equation (

	Remark 3. Example 2.8 could be generalized by considering a measure, say µ, over R, to handle the
	potential contribution of discrete-time events to the pointwise connectivity:
	ż
	Ψ Z,µ pA | t, s, xq " Zps, xq	1 pΦpv,s,xqPAq Zpv, Φpv, s, xqqdµpvq.	(9)
	∆ts
	Remark 4.

located at x at time s) is lower than a threshold h when it is located at Φpv, s, xq at time v; (ii) Z is a function of the local rainfall intensity at ps, xq. Remark 2.

Table 1 :

 1 Network indices (Diameter, density, transitivity, shortest path (mean and standard deviation), small worldness, scale-free property, degree correlation) calculated from the networks covering the Mediterranean region and estimated in three temporal contexts: the entire period 2011-2017, yearly time periods from 2011 to 2017 and monthly time periods.

			Dens Trans S P (mean)	S P (sd)	S W S F D C
	2011-2017 3.12ˆ10 ´3 2.96ˆ10 ´4	0.74	8.20ˆ10 ´4 2.19ˆ10 ´4 18242 11.6	0.31
	2011	0.03	3.0ˆ10 ´4	0.67	6.17ˆ10 ´3 1.92ˆ10 ´3	2222 11.5 -0.04
	2012	0.02	2.9ˆ10 ´4	0.67	6.09ˆ10 ´3 1.87ˆ10 ´3	2236 10.4 -0.06
	2013	0.01	3.0ˆ10 ´4	0.66	6.06ˆ10 ´3 1.86ˆ10 ´3	2228 13.9	0.20
	2014	0.08	2.9ˆ10 ´4	0.65	6.52ˆ10 ´3 2.11ˆ10 ´3	2044	9.0	0.30
	2015	0.06	3.0ˆ10 ´4	0.66	6.22ˆ10 ´3 1.93ˆ10 ´3	2147	8.2	0.26
	2016	0.55	3.0ˆ10 ´4	0.66	6.32ˆ10 ´3 1.97ˆ10 ´3	2110	8.8	0.25
	2017	0.12	3.0ˆ10 ´4	0.65	6.08ˆ10 ´3 1.83ˆ10 ´3	2186	9.2	0.12
	January	0.23	3.1ˆ10 ´4	0.63	1.15ˆ10 ´2 3.81ˆ10 ´3	1118 14.9	0.14
	February	0.67	3.0ˆ10 ´4	0.63	1.15ˆ10 ´2 3.92ˆ10 ´3	1106 14.9	0.24
	March	0.29	3.0ˆ10 ´4	0.62	1.15ˆ10 ´2 3.87ˆ10 ´3	1104 11.9	0.26
	April	0.70	3.1ˆ10 ´4	0.64	1.13ˆ10 ´2 3.91ˆ10 ´3	1137 15.4	0.18
	May	1.02	3.2ˆ10 ´4	0.63	1.20ˆ10 ´2 4.43ˆ10 ´3	1064 12.4	0.21
	June	1.14	3.2ˆ10 ´4	0.61	1.19ˆ10 ´2 4.18ˆ10 ´3	1041 12.9	0.15
	July	1.41	3.1ˆ10 ´4	0.60	1.19ˆ10 ´2 4.01ˆ10 ´3	1017 28.2 -0.07
	August	1.87	2.9ˆ10 ´4	0.59	1.22ˆ10 ´2 4.14ˆ10 ´3	987 11.9 -0.02
	September	1.51	2.9ˆ10 ´4	0.60	1.23ˆ10 ´2 4.18ˆ10 ´3	994 12.8	0.13
	October	0.12	2.8ˆ10 ´4	0.62	1.26ˆ10 ´2 4.46ˆ10 ´3	998 14.9	0.12
	November	0.10	2.7ˆ10 ´4	0.62	1.27ˆ10 ´2 4.90ˆ10 ´3	997 14.3	0.41
	December	1.07	3.0ˆ10 ´4	0.61	1.16ˆ10 ´2 3.89ˆ10 ´3	1079 14.1	0.32

Table 2 :

 2 Network indices (Diameter, density, transitivity, shortest path (mean and standard deviation), small worldness, scale-free property, degree correlation) calculated from the networks covering PACA and estimated in three temporal contexts: the entire period 2011-2017, yearly time periods from 2011 to 2017 and monthly time periods.

Acknowledgments

This research was funded by the SPREE project from the French National Research Agency (grant n 0 ANR-17-CE32-0004-01) and the PHYTOSENTINEL project (grant n 0 IB-2019-SPE). The authors thank Loïc Houde for his technical assistance in the calculation of trajectories with HYSPLIT.

Discussion

We presented a framework for estimating and characterizing spatial and spatio-temporal networks generated by trajectory data. The development of this framework was motivated by the study of networks resulting from the movement of air masses sampled over long time periods and large spatial scales. Thus, in the application, we investigated the tropospheric connectivities across the Mediterranean basin and the French region PACA, and their variations through years and months. Our approach could be applied to diverse phenomena, from which trajectories can be observed. For instance, one could estimate networks generated by the movement of animals on the landscape scale based on animal trajectories observed with GPS devices [START_REF] Bastille-Rousseau | Applying network theory to animal movements to identify properties of landscape space use[END_REF]. This would allow the characterization of connectivity between different landscape components. Sampled trajectories of humans, sampled transports of specific goods (such as plant material) and sampled trajectories of knowledge in social communities (that cannot be exhaustively observed) could also be used to estimate networks in other applied settings.

In Section 2.3, we proposed diverse measures of connectivity with different underlying (physical or biological) interpretations. Thus, the analyst can adapt the connectivity measure to the mechanistic processes he investigates. In the application section, we only used the contact-based connectivity. Comparisons of contact-based, length-based and duration-based connectivities, not shown in this manuscript, led to little variations for the two case studies considered in this article. However, the use of covariates such as local rainfall and air-mass altitude for defining connectivities, as proposed in Section 2.3, is expected to potentially impact the inferred networks and deserves to be explored. This perspective would be particularly relevant in the context of aerobiology: e.g., the airborne transport of organic particles, such as bacteria and fungal spores, can be influenced by rainfall favoring the deposition of these particles [START_REF] Morris | Mapping rainfall feedback to reveal the potential sensitivity of precipitation to biological aerosols[END_REF].

In statistics, we are not only interested in point estimation, but also in the assessment of estimation uncertainties. In this paper, we however, focused on connectivity estimation, even if quantifying the estimation variance could have been useful for more rigorously investigating temporal variation in connectivities. Formally, the connectivity measures that we defined are integrals. Hence, results on integral numerical approximations (e.g., midpoint, trapezoidal or Monte Carlo integration) can be exploited to assess errors or variances of the connectivity estimates [START_REF] Davis | Methods of numerical integration[END_REF][START_REF] Caflisch | Monte carlo and quasi-monte carlo methods[END_REF][START_REF] Geweke | Monte carlo simulation and numerical integration[END_REF]. However, for this assessment, one should ideally take into account dependencies between connectivity estimates for different pairs of nodes, which is not trivial. Further in-depth methodological developments are required to tackle this issue.

To more finely estimate connectivity, and its uncertainty, one could also take into account, if relevant, the uncertainty about the trajectories themselves. For example, when observed trajectories are smoothed versions of actual trajectories (as it is likely the case for air-mass trajectories calculated with HYSPLIT)

or when the trajectories are partially observed and rather erratic, (i) a probabilistic model grounded on, for instance, a stochastic differential equation, could be used to reconstruct probable trajectories and (ii) the connectivity would be estimated from these reconstructed trajectories. Obviously, step (ii) should