Zhouchuan Huang
email: huangz@nju.edu.cn

Sébastien Chevrot

Mantle dynamics in the SE Tibetan Plateau revealed by teleseismic shearwave splitting analysis

Keywords: Shear-wave splitting, Seismic tomography, ChinArray, Two-layer anisotropy model, Tibetan Plateau, Asthenospheric extrusion

We measure frequency-dependent shear wave splitting parameters in the SE Tibetan Plateau.

2. Two-layer anisotropy better explains the observations at many stations.

3. The asthenosphere beneath the plateau is extruded through the SE Tibetan Plateau.

Introduction

The Tibetan Plateau results from the collision between the Indian and Eurasian Plates which started ~ 50 Ma (Figure 1) [START_REF] Yin | Geologic evolution of the Himalayan-Tibetan orogen[END_REF][START_REF] Tapponnier | Oblique stepwise rise and growth of the Tibet plateau[END_REF][START_REF] Royden | The geological evolution of the Tibetan plateau[END_REF]. It is a major Cenozoic structure that affects the environment and climate of whole Asia [START_REF] Yin | Geologic evolution of the Himalayan-Tibetan orogen[END_REF], which makes it a major focus in geosciences. Many different models were proposed to explain the growth of the Tibetan Plateau. One class of models suggests that the plateau is caused by the deformation of the whole lithosphere, produced either by plate subduction [START_REF] Tapponnier | Oblique stepwise rise and growth of the Tibet plateau[END_REF][START_REF] Kind | Seismic images of crust and upper mantle beneath Tibet: Evidence for Eurasian plate subduction[END_REF] or distributed lithospheric deformation [START_REF] England | Finite strain calculations of continental deformation: 2. Comparison with the India!Asia Collision Zone[END_REF]. Another kind of models emphasizes crustal deformations, in which the plateau uplift is produced by thickened upper crust (due to pure shear deformations) [START_REF] Zhang | Continuous deformation of the Tibetan Plateau from global positioning system data[END_REF] or lower crust (due to crustal flow) [START_REF] Royden | Surface deformation and lower crustal flow in eastern Tibet[END_REF][START_REF] Clark | Topographic ooze: Building the eastern margin of Tibet by lower crustal flow[END_REF].

The SE Tibetan Plateau is key to understanding the tectonic evolution of the high plateau. It is bounded by the stable Sichuan Basin to the northeast and the eastern Himalayan syntaxis (EHS) to the southwest (Figure 1). GPS measurements indicate that the upper crust is rotating clockwise around the EHS [START_REF] Zheng | Crustal Deformation in the India!Eurasia Collision Zone From 25 Years of GPS Measurements[END_REF]. The derived crustal strains are generally localized near major faults such as the Lijiang-Xiaojinhe fault at the plateau margin and the large sinistral Xianshuihe-Xiaojiang fault [START_REF] Kreemer | A geodetic plate motion and Global Strain Rate Model[END_REF][START_REF] Zheng | Crustal Deformation in the India!Eurasia Collision Zone From 25 Years of GPS Measurements[END_REF].

However, focal mechanisms indicate that the stress field in the upper crust is controlled by the gravity potential [START_REF] Xu | Crustal stress field in Yunnan: implication for crust-mantle coupling[END_REF][START_REF] Xu | Sharp lateral Moho variations across the SE Tibetan margin and their implications for plateau growth[END_REF].

Many seismological studies have been conducted to reveal the deep architecture of the SE Tibetan Plateau. The crustal thickness decreases from more than 60 km in eastern Tibet to 30 km in the surrounding Yangtze and Indochina blocks to the southeast [START_REF] Wang | Crustal thickness and Poisson's ratio in southwest China based on data from dense seismic arrays[END_REF][START_REF] Xu | Sharp lateral Moho variations across the SE Tibetan margin and their implications for plateau growth[END_REF]. Seismic velocities in the mid-lower crust beneath the plateau are generally lower than those in the surrounding regions [START_REF] Yao | Heterogeneity and anisotropy of the lithosphere of SE Tibet from surface wave array tomography[END_REF][START_REF] Liu | Eastward expansion of the Tibetan Plateau by crustal flow and strain partitioning across faults[END_REF]Bao et al., 2015b]. Recent high-resolution tomography revealed two belts of low seismic velocity, high attenuation and conductivity in the lower crust, which were interpreted as channel flows in the deep crust of the SE Tibetan Plateau [Bao et al., 2015a]. Similar channel flows may also exist in the asthenosphere. Regional and teleseismic tomography imaged a low velocity zone in the upper mantle surrounding the southwestern margin of the stable, high velocity root of the Yangtze craton [Huang et al., 2015a;2015c;[START_REF] Lei | Teleseismic P-wave tomography and mantle dynamics beneath Eastern Tibet[END_REF]. Mantle processes may also be affected by the subduction of the Indian Plate, which was imaged as a high velocity anomaly that extends from the uppermost mantle to the mantle transition zone [START_REF] Li | Subduction of the Indian lithosphere beneath the Tibetan Plateau and Burma[END_REF][START_REF] Wei | East Asia mantle tomography: New insight into plate subduction and intraplate volcanism[END_REF][START_REF] Xu | Insight into the subducted Indian slab and origin of the Tengchong volcano in SE Tibet from receiver function analysis[END_REF].

Flow in the mantle produces a lattice preferred orientation of olivines and pyroxenes, and thus seismic anisotropy, which may be captured by seismic waves that travel through it [START_REF] Silver | Shear-Wave Splitting and Subcontinental Mantle Deformation[END_REF][START_REF] Mainprice | Seismic Anisotropy of the Deep Earth from a Mineral and Rock Physics Perspective[END_REF][START_REF] Maupin | Theory and Observations -Seismic Anisotropy[END_REF]. Shear waves split into two quasi shear waves when they travel through an anisotropic medium [START_REF] Silver | Shear-Wave Splitting and Subcontinental Mantle Deformation[END_REF]. These shear waves are orthogonally polarized and propagate with different velocities, a physical phenomenon called birefringence or shear wave splitting. Two parameters are used to describe this shear-wave splitting quantitively. The polarization of the fast wave (FPD, ) indicates the orientation of the aligned fabrics in the underlying medium. The delay time ( ) is the time-difference accumulated by the two quasi shear waves. It depends on the strength of anisotropy and on the thickness of the anisotropic layer. Crustal anisotropy beneath the SE Tibetan Plateau has been investigated intensively [START_REF] Yao | Heterogeneity and anisotropy of the lithosphere of SE Tibet from surface wave array tomography[END_REF][START_REF] Shi | Shear-wave splitting beneath Yunnan area of Southwest China[END_REF][START_REF] Sun | Crustal structure and deformation of the SE Tibetan plateau revealed by receiver function data[END_REF][START_REF] Sun | Effects of olivine fabric, melt-rock reaction, and hydration on the seismic properties of peridotites: Insight from the Luobusha ophiolite in the Tibetan Plateau[END_REF][START_REF] Cai | Crustal anisotropy and deformation of the southeastern margin of the Tibetan Plateau revealed by Pms splitting[END_REF][START_REF] Huang | P Wave Anisotropic Tomography of the SE Tibetan Plateau: Evidence for the Crustal and Upper-Mantle Deformations[END_REF][START_REF] Han | Layered crustal anisotropy and deformation in the SE Tibetan plateau revealed by Markov-Chain-Monte-Carlo inversion of receiver functions[END_REF]. The orientation of microcracks in the upper crust, revealed mainly by local shear wave splitting and surfacewave studies, are usually parallel to the strike of faults [START_REF] Yao | Heterogeneity and anisotropy of the lithosphere of SE Tibet from surface wave array tomography[END_REF][START_REF] Shi | Shear-wave splitting beneath Yunnan area of Southwest China[END_REF][START_REF] Huang | P Wave Anisotropic Tomography of the SE Tibetan Plateau: Evidence for the Crustal and Upper-Mantle Deformations[END_REF][START_REF] Han | Layered crustal anisotropy and deformation in the SE Tibetan plateau revealed by Markov-Chain-Monte-Carlo inversion of receiver functions[END_REF]. The anisotropy in the lower crust is different. Splitting measurements performed on Pms phases indicate large delay times [START_REF] Sun | Crustal structure and deformation of the SE Tibetan plateau revealed by receiver function data[END_REF][START_REF] Cai | Crustal anisotropy and deformation of the southeastern margin of the Tibetan Plateau revealed by Pms splitting[END_REF]. Since upper crust accounts for up to 0.1 s delay times, the Pms delay times (~0.3 s) mostly reflect anisotropy in the lower crust [START_REF] Han | Layered crustal anisotropy and deformation in the SE Tibetan plateau revealed by Markov-Chain-Monte-Carlo inversion of receiver functions[END_REF]. The Pms splitting FPDs beneath the high plateau are consistent with the model of deep crustal flow from the plateau to its southeast margin [START_REF] Sun | Crustal structure and deformation of the SE Tibetan plateau revealed by receiver function data[END_REF]. In contrast, recent high-resolution P wave anisotropic tomography found that the fast orientation is parallel to the plateau margin as well as the contour lines of the relief and Moho depths [START_REF] Huang | P Wave Anisotropic Tomography of the SE Tibetan Plateau: Evidence for the Crustal and Upper-Mantle Deformations[END_REF][START_REF] Han | Layered crustal anisotropy and deformation in the SE Tibetan plateau revealed by Markov-Chain-Monte-Carlo inversion of receiver functions[END_REF]Xu et al., 2020]. These results would indicate that crustal flow has been blocked by the strong Yangtze and Indochina crustal blocks [START_REF] Huang | P Wave Anisotropic Tomography of the SE Tibetan Plateau: Evidence for the Crustal and Upper-Mantle Deformations[END_REF][START_REF] Xu | Sharp lateral Moho variations across the SE Tibetan margin and their implications for plateau growth[END_REF].

High-resolution seismic anisotropy in the upper mantle can be revealed by teleseismic shear wave (e.g., SKS, PKS, SKKS, hereafter XKS for short) splitting [START_REF] Lev | Seismic anisotropy in Eastern Tibet from shear wave splitting reveals changes in lithospheric deformation[END_REF][START_REF] Sol | Geodynamics of the southeastern Tibetan Plateau from seismic anisotropy and geodesy[END_REF][START_REF] Wang | Evidence for mechanically coupled lithosphere in central Asia and resulting implications[END_REF][START_REF] Huang | Three-dimensional P wave azimuthal anisotropy in the lithosphere beneath China[END_REF][START_REF] Chang | Upper mantle anisotropy beneath the southern segment of North-South tectonic belt[END_REF]Huang et al., 2015b;φ δt Tiwari et al., 2017;[START_REF] Kong | Azimuthal anisotropy and mantle flow underneath the southeastern Tibetan Plateau and northern Indochina Peninsula revealed by shear wave splitting analyses[END_REF]. The first-order pattern of XKS splitting is the clockwise rotation in the eastern Tibetan Plateau and the abrupt change of FPDs across ~26°N, from N-S in the north to E-W in the south of the SE Tibetan Plateau. This transition may be explained by a transition in the mode of deformation, from simple shear beneath the high plateau to pure shear beneath the plateau margin [START_REF] Wang | Evidence for mechanically coupled lithosphere in central Asia and resulting implications[END_REF]. However, the thin lithosphere beneath south Yunnan (< 80 km) cannot produce the observed delay times (~1.0-1.5 s), therefore anisotropy in the asthenosphere is also required [Huang et al., 2015b].

The change of the FPDs across 26°N may also be explained by the asthenospheric extrusion that surrounds the southwestern margin of the stable Yangtze craton [Huang et al., 2015b;2015c;[START_REF] Kong | Azimuthal anisotropy and mantle flow underneath the southeastern Tibetan Plateau and northern Indochina Peninsula revealed by shear wave splitting analyses[END_REF]. But all these studies are based on the assumption of one-layer anisotropy with horizontal symmetry axis. The potential complex anisotropic structures such as multi-layer anisotropy as discussed in the Tibetan Plateau are usually ignored. Another potential complex anisotropy is the frequency-dependent shear-wave splitting identified in subduction zones [START_REF] Huang | Frequency-dependent shear-wave splitting and multilayer anisotropy in northeast Japan[END_REF][START_REF] Marson!pidgeon | Frequency!dependent anisot-ropy in Wellington, New Zealand[END_REF][START_REF] Wirth | Frequency!dependent shear wave split-ting beneath the Japan and Izu!Bonin subduction zones[END_REF], which may be attributed to heterogenous anisotropic structures at various scales. [START_REF] Huang | Insight into NE Tibetan Plateau expansion from crustal and upper mantle anisotropy revealed by shear-wave splitting[END_REF] made similar observations in the northeastern Tibetan Plateau.

Poor vertical resolution of XKS splitting and coarse station distribution in this region prevent us from constraining the depth distribution of deformation in the upper mantle. In this study, we measure shear-wave splitting parameters on the XKS waves recorded by the ChinArray experiment deployed in the SE Tibetan Plateau. The measurements are analyzed in different frequency bands and fit by two-layer anisotropy models to detect possible complex anisotropic structures, especially the depth-dependent distribution of shear wave anisotropy in the upper mantle.

Data and method

Data

We collected waveforms recorded by 350 broad band stations of the ChinArray project deployed in the SE Tibetan Plateau during 2011-2013 (Figure 2a). The stations are equipped with a Guralp CMG-3EPC three-component broadband seismometer and a Reftec-130 digitizer (with a sampling rate of 100 samples per second). We extracted records of events with magnitude larger than 5.8 and epicentral distances of 85° -150°, from which we could obtain 173 good events (Figure 2b). A large number of events occurred in Tonga and New Zealand subduction zones. The SKS phases are clear for these earthquakes at distances of 85°-110°, especially for the deep earthquakes in the mantle transition zone. Other groups of events located in North and Central America provided SKKS and PKS phases that could also be exploited at distances of 120-135°.

Shear-wave splitting measurements with minimum-transverse-energy method

In this study, we first used the SplitLab toolbox [START_REF] Wüstefeld | SplitLab: A shear-wave splitting environment in Matlab[END_REF] and the transverse-component minimization method [START_REF] Silver | Shear-Wave Splitting and Subcontinental Mantle Deformation[END_REF] to measure shear wave splitting parameters. The waveforms are filtered by a Butterworth broadband filter of 2-8 s (HF) and 8-20 s (LF), respectively. A grid-search approach is applied to find the optimal splitting parameters that minimize the energy on the corrected transverse component, which was constructed by rotating (-90° to 90° at 1° step) and time-shifting (0 s to 4 s at 0.02 s step) the Q-T components in a ray-coordinate based L-Q-T coordinate system. The uncertainties are estimated by the 95% confidence region of the F-test.

We follow the semi-automatic procedure of [START_REF] Huang | Insight into NE Tibetan Plateau expansion from crustal and upper mantle anisotropy revealed by shear-wave splitting[END_REF] to select the time windows for shear wave splitting analysis (Figure 3). We determine the beginning of the time window from the STA/LTA ratio (short-time to long-time average ratio) on the radial component around the theoretical arrival time of the XKS phases. The end of the time window changes iteratively from 1 to 2.5 times (with a step of 1/10 of) the dominant period after the beginning. The iterations stop if the ending falls into the window of the direct S wave or if the energy on the transverse component increases (see [START_REF] Huang | Insight into NE Tibetan Plateau expansion from crustal and upper mantle anisotropy revealed by shear-wave splitting[END_REF] for details). Finally, after a visual check we rejected poor measurements (e.g., noisy data, poor correction). This semi-automatic procedure improves the reliability of the results [Huang et al., 2017] and makes it convenient to analyze shear wave splitting in different frequency bands.

Two-layer model analysis

We applied the grid-search algorithm of [START_REF] Silver | The Interpretation of Shear!Wave Splitting Parameters In the Presence of Two Anisotropic Layers[END_REF] to search for the optimal two-layer anisotropy models, including the FPDs and time delays of the lower (

) and upper layers ( ). The synthetic splitting parameters () can thus be calculated [START_REF] Silver | The Interpretation of Shear!Wave Splitting Parameters In the Presence of Two Anisotropic Layers[END_REF]. The misfit between the synthetic and observations parameters is minimized for the optimal models, where ( ) and ( ) are the i-th ( ) observational and synthetic parameters, respectively, and are the corresponding 95% uncertainties. and change between 0º-180º with a step of 1º; and change between 0-2.0 s with a step of 0.05 s in the grid search. The improvement of the two-layer models over single-layer models assuming horizontal axes [START_REF] Fontaine | Upper-mantle flow beneath French Polynesia from shear wave splitting[END_REF][START_REF] Huang | Insight into NE Tibetan Plateau expansion from crustal and upper mantle anisotropy revealed by shear-wave splitting[END_REF] is quantified by

φ L , δt L φ U , δt U = N ∑ i=1 ( φ obs i -φ syn i σ φ i ) 2 + ( δt obs i -δt syn i σ δt i ) 2 φ L φ U δt L δt U
where is the number of observations (i.e., 2*N), and k(=4) is the number of model parameters,

) is the reduction of the misfits between observational and synthetic splitting parameters for the two-layer models ( ) compared with the optimal onelayer model ( ). We further derived the 95% uncertainties of the two-layer models such that the models generating the misfit lower than ( is the inverse of the F-distribution) are considered reliable. Parameter provides a statistic estimate to the improvement of the two-layer model although no accurate uncertainties are derived. For example means the two-layer models explain more than 20% of the azimuthal variation of the splitting parameters [START_REF] Fontaine | Upper-mantle flow beneath French Polynesia from shear wave splitting[END_REF][START_REF] Walker | Shear-wave splitting around the Eifel hotspot: evidence for a mantle upwelling[END_REF]; in this study we used 0.2 as the critical value for favoring a two-layer model over a single-layer model.

2.4.

Results

Frequency-dependent shear wave splitting

The XKS phases are clear in both low (8-20 s) and high (2-8 s) frequencies (Fig. 4).

The signal-to-noise ratios are mostly greater than 10 for the "Good" and "Fair" measurements. The splitting parameters are better constrained in the high frequency than in the low frequency (Fig. 4). While the uncertainties of the FPD and delay times in the high frequency are usually smaller than 20° and 0.6 s, respectively.

r min 2 *(1 + f (4,N d -4)*4/(N d -4)) f R* R* ≥ 0.2 ()() ()
Figure 5 shows the shear wave splitting measurements with the traditional minimumtransverse-energy method [START_REF] Silver | Shear-Wave Splitting and Subcontinental Mantle Deformation[END_REF] at low frequency (8-20 s), including 435 "Good" and 768 "Fair" results among a total of 9053 useful measurements after analyzing a total of 42370 records. The FPDs show a dominant transition from N-S to E-W southward across ~26°N. The delay times are mostly between 0.8 and 2.0 s. Strong splitting delays are generally measured at stations in the southern margin of the Sichuan Basin and in the Youjiang Orogen. These shear wave splitting measurements are close to those reported in previous studies. Some differences can be noted, however, in the Tibetan Plateau (e.g., region

A with NE-SW FVDs), as already noted in several previous studies [START_REF] Wang | Evidence for mechanically coupled lithosphere in central Asia and resulting implications[END_REF][START_REF] Chang | Upper mantle anisotropy beneath the southern segment of North-South tectonic belt[END_REF][START_REF] Bao | Episodic Lithospheric Deformation in Eastern Tibet Inferred From Seismic Anisotropy[END_REF]. We also obtained 1046 "Good" and 2626 "Fair" nulls.

Their back azimuths concentrate around 30° and 120° directions, where most earthquakes occurred.

Figure 6 shows the shear wave splitting measurements at higher frequency (2-8 s), including 1108 "Good" and 1538 "Fair" measurements as well as 57 "Good" and 423 "Fair" nulls. The angular differences between the measured FPDs in the higher and lower frequencies are smaller than 20° (Figure 7). However, the delay times differ significantly; the high frequency (~0.5-1.4 s) measurements are much smaller than those at lower frequency (Figure 7c). Among the 497 splitting parameters measured in both higher (2-8 s) and lower (8-20 s) frequencies, the differences of 155 delay times (over 30%) are larger than 0.5 s.

Frequency-dependent shear-wave splitting has been observed in subduction zones [START_REF] Huang | Frequency-dependent shear-wave splitting and multilayer anisotropy in northeast Japan[END_REF][START_REF] Marson!pidgeon | Frequency!dependent anisot-ropy in Wellington, New Zealand[END_REF][START_REF] Wirth | Frequency!dependent shear wave split-ting beneath the Japan and Izu!Bonin subduction zones[END_REF]. The differences in delay times can be explained by the frequency-dependent sensitivities of XKS waves to the anisotropic structures at different scales [START_REF] Marson!pidgeon | Frequency!dependent anisot-ropy in Wellington, New Zealand[END_REF]]. In the viewpoint of finite frequency, lower frequencies sample a wider region and yield a regional average [START_REF] Eakin | Upper mantle anisotropy beneath Peru from SKS splitting: Constraints on flat slab dynamics and interaction with the Nazca Ridge[END_REF][START_REF] Monteiller | High-resolution imaging of the deep anisotropic structure of the San Andreas Fault system beneath southern California[END_REF]Lin et al., 2014aLin et al., , 2014b]]. Figure 8 shows the sensitivity kernels of the splitting-intensity parameters for an XKS waves at periods of 10 s (a, b) and 5 s (c, d), respectively. At 150 km depth, the diameter of the sensitivity kernel at 10-s is 1.5 times that of at 5-s. In addition, forward modeling showed that the measurements in the higher frequency range are generally biased towards anisotropy at shallow depth [START_REF] Saltzer | How are vertical shear wave splitting measurements affected by variations in the orientation of azi-muthal anisotropy with depth[END_REF].

Our XKS splitting measurements differ from crustal anisotropy characterized from the analysis of splitting on S and Pms phases. Local S wave splitting shows that the upper crustal FPDs are parallel to the fault strikes [START_REF] Shi | Shear-wave splitting beneath Yunnan area of Southwest China[END_REF], while Pms splitting revealed different deep crustal anisotropy due to specific lower crustal deformations during the plateau expansion [START_REF] Sun | Crustal structure and deformation of the SE Tibetan plateau revealed by receiver function data[END_REF][START_REF] Cai | Crustal anisotropy and deformation of the southeastern margin of the Tibetan Plateau revealed by Pms splitting[END_REF]. Crustal anisotropy only generates splitting delay times up to ~0.3 s, which represents only a small part of the XKS delay times (1-1.5 s) at low frequency but can account for 50% of the XKS delay times at high frequency.

Consequently, XKS splitting especially for those in lower frequency is mainly caused by anisotropy in the upper mantle. The strong crustal small-scale heterogeneity may produce only a small coherent contribution to apparent splitting at low frequency.

Two-layer models

We fit the shear wave splitting measurements at every station with two-layer anisotropy models; Figure 9 shows an example. The two-layer models can better explain the azimuthal variations of splitting parameters in the higher frequency band at 28 stations (Figure 10, blue and red rose diagrams). However, the two-layer models only significantly improve the measurements in the lower frequency at 3 stations (Figure 10, cyan and pink rose diagrams); they are similar to the two-layer models for high-frequency measurements at the adjacent stations. As mentioned above, the low-frequency XKS waves have large sensitivity kernels and they capture the longer wavelength anisotropic structures. Thus lateral structural variations may destroy the possible azimuthal variation of the splitting parameters produced by a two-layer model.

In these two-layer models, the FPDs of the lower layer are generally NW-SE.

Nevertheless, at station 53033 to the northeast of the Tengchong volcano, the lower-layer

FPDs are NEE-SWW to nearly E-W. [START_REF] Kong | Azimuthal anisotropy and mantle flow underneath the southeastern Tibetan Plateau and northern Indochina Peninsula revealed by shear wave splitting analyses[END_REF] analyzed the shear splitting at permanent stations in the SE Tibetan Plateau. They obtained the similar E-W FPDs in the lower layer at stations in the northern Indochina and the Red-River fault zone (Figure 10).

The FPDs of the upper layer show strong lateral variations (Figure 10). They are NW-SE to N-S and parallel to the strikes of faults in the Tibetan Plateau. In the southern margin of the Sichuan Basin, the upper-layer and lower-layer FPDs are not well discriminated actually, both NW-SE. At station 53231, the upper-layer FPDs rotate to be nearly N-S, which are consistent to the fault strikes. In the northern Indochina and Red-River fault zone, the FPDs are nearly N-S. In the Youjiang Orogen, the FPDs are mostly E-W.

3.3.

Discussion

Lithospheric anisotropy

The SE Tibetan Plateau consists of many blocks with different sizes and vertical extents. To first order, the tectonics of the region is dominated by the interactions between the stable Yangtze craton, the Indochina blocks, and a series of sub-blocks belonging to the Tibetan Plateau (Figure 2a) [START_REF] Royden | The geological evolution of the Tibetan plateau[END_REF]. The Youjiang Orogen is amalgamated with the Yangtze craton along the southern margin and collided with the Indochina block along the Red-River faults zone during the closure of the Paleo-Tethyan ocean in the Mesozoic [START_REF] Yang | Detrital record of Indosinian mountain building in SW China: Provenance of the Middle Triassic turbidites in the Youjiang Basin[END_REF][START_REF] Faure | The South China block-Indochina collision: Where, when, and how[END_REF][START_REF] Qiu | Mesozoic geology of southwestern China: Indosinian foreland overthrusting and subsequent deformation[END_REF]. The southeast Chuandian block (among Lijiang-Xiaojinhe, Red-River, and Xiaojiang faults) is the southwestern part of the Yangtze craton but it has been involved in the active tectonics of the SE Tibetan Plateau.

The present study only covers the SE margin of the Tibetan Plateau. The shear-wave splitting measurements show fault-parallel FPDs in general (Figures 567), in particular in the upper-layer of two-layer models (Figure 10). Thus they mainly arise from the lithospheric deformations caused by the northeast subduction/collision of the eastern Himalayan syntax.

Along the southeastern plateau margin, there are many NE-SW FPDs, which are parallel to the plateau margin (A and B in Figures 567). These results obviously reflect the anisotropy in the shallow part, as can be inferred from our synthetic tests. They are therefore caused by the lithospheric deformations due to the interaction between the Tibet Plateau and the Yangtze craton [e.g., [START_REF] Bao | Episodic Lithospheric Deformation in Eastern Tibet Inferred From Seismic Anisotropy[END_REF].

The lithospheric anisotropy of the Yangtze craton is relatively simple. The lithospheric thickness in the core of the Yangtze craton beneath the Sichuan Basin is ~200 km [START_REF] Pasyanos | LITHO1.0: An updated crust and lithospheric model of the Earth[END_REF]. The NW-SE FPDs spread in the whole lithosphere (Figures 567), mostly representing the fossil anisotropy frozen in the lithosphere. Similar processes may explain the E-W FPDs along the southern margin of the Yangtze carton where the Youjiang

Orogen is amalgamated to form the South China block in the Triassic [START_REF] Yang | Detrital record of Indosinian mountain building in SW China: Provenance of the Middle Triassic turbidites in the Youjiang Basin[END_REF][START_REF] Faure | The South China block-Indochina collision: Where, when, and how[END_REF][START_REF] Qiu | Mesozoic geology of southwestern China: Indosinian foreland overthrusting and subsequent deformation[END_REF]. The E-W FPDs are generally parallel to the boundary between the Youjiang Orogen and Yangtze craton and likely developed during the block amalgamation.

The Red-River fault zone along the southwestern margin of the South China block (including Yangtze craton and Youjiang Orogen) was previously part of the Mesozoic suture zone between the South China and Indochina blocks [START_REF] Faure | The South China block-Indochina collision: Where, when, and how[END_REF]. It was then reactivated as a strike-slip fault during the Indo-Asian collision and the simultaneous plateau growth [START_REF] Tapponnier | Propagating Extrusion Tectonics in Asia -New Insights From Simple Experiments with Plasticine[END_REF][START_REF] Tapponnier | Oblique stepwise rise and growth of the Tibet plateau[END_REF]. The fault-parallel FPDs (i.e., NW-SE) especially of the upper layer (Figure 10) well reflect the aligned mineral fabrics along the tectonic boundary. The weak lithospheric anisotropy beneath the northern Indochina block may result from the small thickness of the lithosphere ~60-80 km [START_REF] Pasyanos | LITHO1.0: An updated crust and lithospheric model of the Earth[END_REF][START_REF] Hu | Seismic upper mantle discontinuities beneath Southeast Tibet and geodynamic implications[END_REF].

Asthenospheric extrusion/flow

The anisotropy of the lower-layer allows us to distinguish two domains. The first domain is the continuous NW-SE FPDs along the southwestern margin of the Yangtze cratonic root (Figure 10). We propose that the NW-SE FPDs represent the asthenospheric extrusion/flow from the Tibetan Plateau. Seismic tomography indicates northward Indian subduction and southward Eurasian subduction beneath the high plateau [START_REF] Li | Subduction of the Indian lithosphere beneath the Tibetan Plateau and Burma[END_REF][START_REF] Wei | East Asia mantle tomography: New insight into plate subduction and intraplate volcanism[END_REF]Huang et al., 2015a]. The asthenospheric space is largely shrunk so that the asthenosphere is thickened vertically while extruded through mantle channelling along the plateau margin. The SE Tibetan Plateau margin is probably the best region to observe the effect of asthenospheric extrusion because the other margins mostly involve stable blocks (e.g., Tarim, Ordos, Sichuan blocks) [START_REF] Huang | High!resolution mantle tomography of China and surrounding regions[END_REF][START_REF] Wei | East Asia mantle tomography: New insight into plate subduction and intraplate volcanism[END_REF]. This is in good agreement with surface wave tomography, which clearly imaged a thickened asthenospheric channel under the SE Tibetan Plateau [START_REF] Zhang | Structural model of the lithosphere-asthenosphere system beneath the Qinghai-Tibet Plateau and its adjacent areas[END_REF]. Another factor leading to the asthenospheric flow is an 80 km step at the bottom of the lithosphere. The lithospheric thickness beneath the Chuandian block and Youjiang Orogen is less than 80 km, in great contrast to >160 km lithosphere beneath the Sichuan Basin (Figure 11a) [START_REF] Pasyanos | Comparison of lower-layer anisotropy with (a) lithospheric thickness[END_REF][START_REF] Hu | Seismic upper mantle discontinuities beneath Southeast Tibet and geodynamic implications[END_REF]. Consequently, the asthenospheric material may be extruded around the southwestern margin of the cratonic root. Seismic tomography revealed much lower velocities around the southwestern cratonic margin (Figure 11b) [Huang et al., 2015c;[START_REF] Huang | lithospheric thickness and P wave velocity anomalies with scale shown at the bottom[END_REF][START_REF] Lei | Is there a big mantle wedge under eastern Tibet?[END_REF], supporting possible asthenospheric extrusion.

The second group of anisotropy is the dominant E-W FPDs extending from northern Indochina block to Youjiang Orogen (Figure 10). The E-W FPDs are mainly caused by the mantle flow related to the eastward subduction of the Indian plate which has been clearly revealed by tomographic studies [START_REF] Wei | East Asia mantle tomography: New insight into plate subduction and intraplate volcanism[END_REF]Huang et al., 2015a;[START_REF] Lei | Is there a big mantle wedge under eastern Tibet?[END_REF]. The northern Indochina and Youjiang blocks are characterized by very thin lithosphere and low velocities in the upper mantle (Figure 11) [START_REF] Pasyanos | LITHO1.0: An updated crust and lithospheric model of the Earth[END_REF]Huang et al., 2015c[START_REF] Huang | P and S Wave Tomography Beneath the SE Tibetan Plateau: Evidence for Lithospheric Delamination[END_REF],

indicating that they are heavily affected by the Indian subduction and subsequent rollback.

The regional mantle-wedge flow may be responsible for the E-W mineral fabrics and anisotropy in the upper mantle [Huang et al., 2015b;[START_REF] Kong | Azimuthal anisotropy and mantle flow underneath the southeastern Tibetan Plateau and northern Indochina Peninsula revealed by shear wave splitting analyses[END_REF]. Geophysical and geochemistry results also suggest lithospheric delamination and asthenospheric upwelling in northern Indochina [START_REF] Chung | Intraplate extension prior to continental extrusion along the Ailao Shan-Red River shear zone[END_REF][START_REF] Lu | Geochemical, Sr-Nd-Pb, and Zircon Hf-O Isotopic Compositions of Eocene-Oligocene Shoshonitic and Potassic Adakite-like Felsic Intrusions in Western Yunnan, SW China: Petrogenesis and Tectonic Implications[END_REF][START_REF] Chen | Delamination of lithospheric mantle evidenced by Cenozoic potassic rocks in Yunnan, SW China: A contribution to uplift of the Eastern Tibetan Plateau[END_REF][START_REF] Wang | Intraplate extension of the Indochina plate deduced from 26 to 24 Ma A-type granites and tectonic implications[END_REF].

The upwelling may in fact break the horizontal mineral alignment resulting in weaker anisotropy. These predictions are consistent with the observed smaller XKS splitting delay times in this region (Figure 567).

The two domains of deep asthenospheric anisotropy merge in the eastern Youjiang Orogen (Figure 11) while eastward flow is observed beneath eastern China. Regional seismic tomography together with XKS splitting in Southeast Asia imply that the mantle flow may reach the northwestern margin of the South China Sea [START_REF] Wei | East Asia mantle tomography: New insight into plate subduction and intraplate volcanism[END_REF]Huang et al., 2015a;2015b]. This regional mantle flow may contribute to the volcanism in Hainan and southern

Vietnam and possibly affect the evolution of the South China Sea (see [Huang et al., 2015a] for details).

Conclusions

1.

We characterized shear wave splitting in the SE Tibetan Plateau 6. The same as Figure 5 for at higher frequency (2-8 s).

  using XKS waves recorded by the dense ChinArray experiment. The measurements show a dominant change of the FPDs across 26°N, i.e., from N-S in the north to E-W in the south, in good agreement with previous studies. We measured the shear-wave splitting in two different frequency bands (2-8 s and 8-20 s) with the Silver and Chan technique, which give similar FPDs but the delay times measured at 2-8 s are smaller than those at 8-20 s. Moreover, the azimuthal variations of the FPDs at many stations can be better explained by two-layer anisotropy models. , the FPDs in the different blocks show different orientations: N-S to NW-SE in the Tibetan Plateau, NW-SE in the Yangtze craton (Sichuan Basin), and E-W in northern Indochina. The distribution of seismic anisotropy seems to be controlled by the main tectonic boundaries, with NE-SW FPDs along the Tibet-Yangtze boundary, NW-SE FPDs along the Indochina-Yangtze boundary, and E-W FPDs along the Yangtze-Youjiang orogens. These results indicate that lithospheric anisotropies are mainly produced by past and present deformations in the lithosphere. In the asthenosphere, the shear wave anisotropy suggests two channels of asthenospheric flow. One channel extends from the Tibetan Plateau to eastern Youjiang, characterized as NW-SE FPDs, indicating southeastern asthenospheric extrusion along the southwestern margin of the Yangtze cratonic root at depth. Another channel extends from northern Indochina to Youjiang Orogen eastward, characterized as E-W FPDs, resulting from the mantle flow driven by the eastern subduction of the Indian plate. (iteratively moved in the dark grey region) of the picked SKS window. (c) Shear wave splitting measurements in 7 iterations. White dots and short bars denote the splitting parameters and their uncertainties, respectively, in each iteration. Gray dots show the estimated polarization direction. Blue symbols denote the best splitting parameters among all the 7 iterations. (d) Initial (dashed) and corrected (solid) particle motions before and after shear wave splitting corrections. (e) Corrected radial (dashed) and transverse (solid) components. (f) The distribution of the energy on the transverse component in grid search. The shaded region in bold line denote the uncertainties defined by 95% confidence in F-test.4. Comparison of dominant periods (a, e), signal-to-noise ratio (b, f), uncertainties of the FPDs (c, g) and delay times (d, h) of the "Good" and "Fair" measurements in the lower (upper panel) and higher (lower panel) frequency measurements. 5. Shear wave splitting measurements obtained using the transverse-component minimization method at lower frequency(8-20 s) [START_REF] Silver | Shear-Wave Splitting and Subcontinental Mantle Deformation[END_REF]. (a) Red and blue bars denote the good and fair non-null measurements. Their orientations and lengths denote the fast polarization directions (FPD) and delay times (scale shown in bottom right), respectively. Yellow bars denote the previous results [Chang et al., 2015]. Labels A and B mark two regions with NE-SW FPDs along the Tibet-Yangtze boundary. The red triangle denotes the Tengchong volcano. The black lines denote the active faults. The black arrow denotes the absolute plate-motion directions of the eastern Himalayan syntaxis (EHS) in the NNR-MORVEL56 model [Argus et al., 2011]. (b) Histograms of the delay times for the good (red) and fair (blue) measurements. (c) Histograms of the back-azimuths for the good (red) and fair (blue) nulls.

Acknowledgements. The waveform data used in this study were provided by the China Seismic Array Data Management Center at the Institute of Geophysics, China Earthquake Administration (ChinArray DMC, doi:10.12001/ChinArray.Data, http:// www.chinarraydmc.cn/). We thank Prof. M. Savage (the Editor), the Associate Editor, Prof.

Tuna Eken, and two anonymous reviewers for their constructive comments and suggestions.

This work was supported by the National Natural Science Foundations of China (41674044) and by the ChinArray Program (DQJB16A0306). ZH is also supported by the Deng-Feng Scholar Program of Nanjing University. Most figures were made using GMT [Wessel et al., 2019].


Diego Vinasco-Alvarez
email: diego.vinasco-alvarez@univ-lyon2.fr

John Samuel
email: john.samuel@cpe.fr

Sylvie Servigne
email: sylvie.servigne@insa-lyon.fr

Gilles Gesquière
email: gilles.gesquiere@univ-lyon2.fr

TOWARDS LIMITING SEMANTIC DATA LOSS IN 4D URBAN DATA SEMANTIC GRAPH GENERATION

Keywords: Graphs, Networks of Ontologies, 3D Urban Data, Spatio-temporal Data, Data Integration, Data Transformation, CityGML, GeoSPARQL

To enrich urban digital twins and better understand city evolution, the integration of heterogeneous, spatio-temporal data has become a large area of research in the enrichment of 3D and 4D (3D + Time) semantic city models. These models, which can represent the 3D geospatial data of a city and their evolving semantic relations, may require data-driven integration approaches to provide temporal and concurrent views of the urban landscape. However, data integration often requires the transformation or conversion of data into a single shared data format, which can be prone to semantic data loss. To combat this, this paper proposes a modelcentric ontology-based data integration approach towards limiting semantic data loss in 4D semantic urban data transformations to semantic graph formats. By integrating the underlying conceptual models of urban data standards, a unified spatio-temporal data model can be created as a network of ontologies. Transformation tools can use this model to map datasets to interoperable semantic graph formats of 4D city models. This paper will firstly illustrate how this approach facilitates the integration of rich 3D geospatial, spatio-temporal urban data and semantic web standards with a focus on limiting semantic data loss. Secondly, this paper will demonstrate how semantic graphs based on these models can be implemented for spatial and temporal queries toward 4D semantic city model enrichment.

INTRODUCTION

The continuing amelioration of semantic 3D city models has provided powerful tools for comprehending the dynamics of the constantly evolving urban landscape. Data-driven approaches such as the construction of virtual environments such as digitaltwins [START_REF] Batty | Digital twins[END_REF][START_REF] Julin | Characterizing 3D City Modeling Projects: Towards a Harmonized Interoperable System[END_REF][START_REF] Schrotter | The Digital Twin of the City of Zurich for Urban Planning[END_REF] can help visualize and simulate city events and dynamics both real and imaginary, over time [START_REF] Biljecki | Applications of 3D City Models: State of the Art Review[END_REF][START_REF] Chaturvedi | A Requirement Analysis on Extending Semantic 3D City Models for Supporting Time-Dependent Properties[END_REF][START_REF] Jaillot | Delivering timeevolving 3D city models for web visualization[END_REF], Samuel et al., 2020).

Currently, the Open Geospatial Consortium's (OGC) 1CityGML2 open data model is often used to facilitate the storage and exchange of virtual 3D city models. In this model, all city objects can be represented in up to five different, well-defined levels-of-detail (LOD0 to LOD4) with increasing accuracy but structural complexity (e.g., a building is composed of walls, roofs, etc.) [START_REF] Kolbe | Representing and Exchanging 3D City Models with CityGML[END_REF]. It defines the three-dimensional geometry, topology, semantics, and appearance of the most relevant topographic objects in the urban context [START_REF] Gröger | CityGML -Interoperable semantic 3D city models[END_REF].

The CityGML standard is still evolving [START_REF] Kutzner | CityGML 3.0: New Functions Open Up New Applications[END_REF] to represent all the complex data requirements of 4D (3D + Time) semantic city models [START_REF] Chaturvedi | A Requirement Analysis on Extending Semantic 3D City Models for Supporting Time-Dependent Properties[END_REF] and may need enrichment through data integration of other heterogeneous data sources to achieve more complete or detailed views of the urban landscape. To facilitate integration, the use of ontologies has been promoted as they provide flexible, machineprocessable formalizations of data models as semantic graphs [START_REF] Claramunt | Ontologies for geospatial information: Progress and challenges ahead[END_REF][START_REF] Psyllidis | Ontology-Based Data Integration from Heterogeneous Urban Systems: A Knowledge Representation Framework for Smart Cities[END_REF]. Ontologies allow then to create semantic urban graphs based on urban data which may require transformation or conversion of data towards a single common data format. However, transformation approaches are naturally prone to semantic data loss and thus limiting this data loss poses a principal integration challenge.

Another challenge -within the context of multisource temporal urban data integration -is the inconsistent use of identifiers of city objects [START_REF] Chaturvedi | Managing Versions and History Within Semantic 3D City Models for the Next Generation of CityGML[END_REF]. Different snapshots of a city from different sources may use different identifiers for the same city objects between concurrent and successive representations of that city. Recently, graph-based change detection approaches (Jaillot et al., 2020, Nguyen andKolbe, 2020) have been proposed to calculate which city objects are being referenced across different city snapshots, even when there are geometric, or semantic differences between the same object in different temporal representations of the city.

As ontologies and graph formats have proven beneficial towards the semantic enrichment of interoperable multidimensional city models, this paper proposes a model-centric ontology-based data integration approach towards limiting semantic data loss in 4D semantic urban data transformations to semantic graph formats. This contribution consists of several aspects: 1. A proposed 4D semantic urban data model, based on a network of ontologies [START_REF] Euzenat | The category of networks of ontologies[END_REF], incorporating the most recent draft of CityGML 3.0 [START_REF] Kutzner | CityGML 3.0: New Functions Open Up New Applications[END_REF], ISO TC211 geographic information standards3 , and a semantic web spatio-temporal query framework.

2. A data transformation approach for the generation of versioned CityGML 3.0 datasets as interoperable semantic RDF graphs. 3. Implementation of these graphs in a spatio-temporal triple store to facilitate geospatial and temporal queries on realworld city data.

The remainder of this paper is structured as follows: Section 2 will describe related works in multidimensional semantic city models and spatio-temporal urban graph implementations; Section 3 will detail the proposed integration approach; Section 4 is dedicated to the implementation of the approach on a realworld city dataset and to discuss its strengths and limitations; Section 5 will conclude the proposed contribution, discuss future works on this research, and provide resources for reproducing the work.

RELATED WORKS

Although research in geometric and semantic data integration of CityGML with Building Information Model (BIM) standards like Industry Foundation Classes (IFC)4 has had much evolution [START_REF] Deng | Mapping between BIM and 3D GIS in different levels of detail using schema mediation and instance comparison[END_REF], temporal 4D semantic city model integration remains a large area of research. As previously stated, many existing standards lack rich, native temporal support, and generally require one of two types of integration: either a data model is extended to integrate temporal data formats or multiple data models are integrated as a single, spatio-temporal data model. This data integration can also take place at two levels: data model integration at the conceptual, physical, or implementation level and integration at the model instance or dataset level. Taking this into consideration, some existing works propose approaches to facilitate this integration for enriching 3D semantic city models.

4D Semantic Urban Data Model Extension

Extending a 3D data model requires starting with a strong, semantic data model. Since its adoption by the OGC in 2012, CityGML has become a commonly used data model standard for this purpose as it is geospatial, easily extensible through Application Domain Extensions (ADEs), and provides a modular semantic conceptual model for various domains of the urban landscape. In the case of CityGML based 4D city models, the conceptual model is traditionally extended to allow better temporal data support [START_REF] Chaturvedi | Managing Versions and History Within Semantic 3D City Models for the Next Generation of CityGML[END_REF]. Some recent works allow the representation of concurrent and successive versions of city models with graph-based versioning frameworks [START_REF] Samuel | Representation of concurrent points of view of urban changes for city models[END_REF] while others allow for time series and time-dependent data representations (such as internet of things sensor data) [START_REF] Biljecki | CityGML Application Domain Extension (ADE): overview of developments[END_REF].

These extensions are usually created manually and may be time-consuming to generate, verify, and update. Due to the utility demonstrated by these proposed models, the next version of CityGML 3.0 standard will incorporate parts of these works as new Versioning and Dynamizer modules.

Ontology-Based 4D Urban Data Integration

When creating a spatio-temporal data model from two or more sources, often a common data format is used to facilitate the integration. For this purpose, much research and official standards have been made for implementing 3D and 4D geospatial data using semantic web technologies, such as RDF (Resource Description Framework) and OWL (Web Ontology Language) [START_REF] Bonduel | Including widespread geometry formats in semantic graphs using RDF literals[END_REF][START_REF] Brink | Linking spatial data: automated conversion of geo-information models and GML data to RDF[END_REF][START_REF] Lemmens | Towards Linked Data and ontology development for the semantic enrichment of volunteered geo-information[END_REF][START_REF] Métral | Extension and Contextualization for Linked Semantic 3D Geodata[END_REF][START_REF] Nuninger | Linking Theories, Past Practices, and Archaeological Remains of Movement through Ontological Reasoning[END_REF][START_REF] Psyllidis | Ontology-Based Data Integration from Heterogeneous Urban Systems: A Knowledge Representation Framework for Smart Cities[END_REF]. These proposals use machine-readable description logics to create ontologies for semantic modeling and integration.

Furthermore, network of ontology-based integration approaches for 4D semantic data integration have also been proposed [START_REF] Hor | A Semantic Graph Database for BIM-GIS Integrated Information Model for an Intelligent Urban Mobility Web Application[END_REF][START_REF] Tran | Semantic Integration of Raster Data for Earth Observation: An RDF Dataset of Territorial Unit Versions with their Land Cover[END_REF][START_REF] Psyllidis | Ontology-Based Data Integration from Heterogeneous Urban Systems: A Knowledge Representation Framework for Smart Cities[END_REF]. These approaches are modular as each ontology integrated into the network can describe a single domain of information or data model. Additionally, extensions to RDF's query language, SPARQL, such as GeoSPARQL5 and stSPARQL6 can be integrated to support queries of 2D geometries and 2D geometric and topological relations, with stSPARQL also supporting temporal queries on graph data.

Works have also been proposed for the automatic generation of geospatial and spatio-temporal semantic ontologies from urban data. Some of these approaches use mapping transformations of conceptual urban data models as UML (Unified Modeling Language) to OWL [START_REF] Brink | Linking spatial data: automated conversion of geo-information models and GML data to RDF[END_REF][START_REF] De Paepe | Automated uml-based ontology generation in oslo2[END_REF]. These approaches implement their own mappings or propose using the standard ISO 19150-27 transformation rules for geospatial data. Other approaches utilize mappings from physical urban data models as XML Schema to OWL [START_REF] Vinasco-Alvarez | From CityGML to OWL[END_REF][START_REF] Vinasco-Alvarez | Towards a semantic web representation from a 3d geospatial urban data model[END_REF][START_REF] Usmani | Automatic Ontology Generation of BIM and GIS Data[END_REF]. These approaches all incur some amount of semantic loss during transformation. In general, UML to OWL transformations may have more well-defined mappings as both modeling languages use similar concepts like classes and properties. However, even the ISO 19150-2 standard can be ambiguous (OGC, 2017) and XML Schema may provide better definitions of spatial semantics [START_REF] Brink | Linking spatial data: automated conversion of geo-information models and GML data to RDF[END_REF].

Summary

In summary, most of these proposals demonstrate effective approaches towards integrating 4D semantic city data for 4D semantic city model enrichment. Implementing strong data standards in these models and preserving their interoperability is key during integration, which requires minimizing semantic loss during transformation. Currently, the upcoming version of CityGML 3.0 is promising to be one of the most popular data models in this context. This paper proposes a model-centric integration approach structured around CityGML 3.0, in its latest draft, and other spatio-temporal data standards. This approach will incorporate network of ontologies integration and automated ontology generation approaches to limit semantic data loss of the original conceptual data models. It will also demonstrate how this model can be used to generate interoperable 4D semantic graphs.

MODEL-CENTRIC DATA INTEGRATION FOR LIMITING SEMANTIC DATA LOSS

This section will define key concepts and the problem, introduce the proposed approach, and demonstrate its implementation.

Problem Definition

Within 4D semantic city models, the entities concerned are often referred to as city objects. These objects can represent buildings, roads, landmarks, and many other urban entities, and are composed of 3 types of data:

• Geospatial: the object's 2D and/or 3D geometric representation, combined with its terrestrial location, often defined by coordinate reference systems. • Temporal: the period when the object exists, represented by two points in time referencing the objects creation and destruction.

• Semantic: data that can describe, categorize, or provide context to the object and its components or their relationships to other objects.

Integration of geospatial and temporal data are the most straightforward as the values that represent these data are real numbers and timestamps. Their integration may require conversion but is only prone to a loss of numeric precision (e.g. loss of geometric or geospatial precision during the conversion of coordinate reference systems [START_REF] Seeger | Spatial referencing and coordinate systems[END_REF]). On the contrary, integrating semantic urban data -either as data models or data -can be more challenging.

Semantic integration at the conceptual level can be done by combining or mapping concepts from 2 or more conceptual models. For this, each conceptual model must be represented by the same modeling language to facilitate integration. Since different modeling languages have varying amounts of expressivity it is possible that a concept may not be able to be represented in its entirety or must be represented differently after transformation. This implies a loss of interoperability and data based on the new model may also share this loss. In this paper, these results are referred to as semantic loss. Once the conceptual models are represented in a homogeneous modeling language, equivalent concepts can be related to one another to integrate the models. Once this conceptual integration is complete, integrated 4D semantic datasets can be generated based on this model to enrich 4D semantic city models.

The following section will introduce the proposed integration methodology for limiting semantic loss in 4D semantic urban data integration through 4D semantic graphs.

4D Semantic Model Integration

For the integration of different 4D semantic urban conceptual models, the approach uses a network of ontologies as a semantic model. To this end the integrated model will be constructed from several existing conceptual models:

• CityGML 3.0 UML model's most recent draft, in particular the Building and Versioning modules for their urban and spatio-temporal semantics. • GML 3.2 model represented through the ISO 19107 standard and GeoSPARQL's GML 3.2 ontology to support CityGML's geometric and geospatial semantics. • GeoSPARQL and stRDF ontologies to provide a spatiotemporal query framework for semantic graph interaction.

These models were chosen for their synergy as a network of ontology integration approach. Most of these models have can be implemented to create OWL ontologies from geospatial UML models. However, implementing this tool for this purpose has been tested by the OGC according to the Testbed 12 ShapeChange Engineering Report (OGC, 2017) and several ambiguities in the ISO 19150-2 standard have been identified which affect the CityGML 3.0 UML model to OWL transformation. The most notable of these being the ambiguous mapping of unions to OWL, which affects CityGML elements such as the union core:cityModelMember. This union groups several important properties of CityGML 3.0 city models, that reference city objects, versions, and version transitions of the model.

To overcome this ambiguity, the Testbed 12 report proposes 3 ShapeChange configurations, of which the "flattening" approach is used for transforming unions:

For each UML property or attribute A -that has a union as its value type -a new property is created for each property option B of the union. The name of the new property is constructed by merging the names of A and B with a union separator as follows:

"[P roperty A][Separator][P roperty B]" (1)
The results of this mapping are shown in figure 2. While this ap- It is also important that the relationships between UML classes and their properties are preserved in the generated ontology. This is ensured by configuring ShapeChange with three mapping transformation rules:

1. Universal quantifications are used in the definition of classes to refer to any properties they may have. 2. Local naming conventions are used to distinguish between different properties of the same name. 3. The domain and range of properties are always declared.

This results in bidirectional references between classes and properties which is useful for automatic semantic graph generation and helps avoid inconsistencies when validating the ontology.

Finally, ShapeChange is configured to map any geometric references to their equivalent in the GeoSPARQL and ISO standard ontologies. Any other configurations follow the OGC Testbed 12 recommendations, enabling ISO 19150-2 rules whenever possible. The final output of the transformation results in 17 linked ontologies, one for each module of the CityGML 3.0 Conceptual Model. Figure 3 exemplifies a mapping between the UML representation of a version and its representation in OWL 2. Next, links are made manually between the ISO 19107 and GeoSPARQL ontologies. Once validated, this is considered as a 4D semantic urban data network of ontologies.

Graph Generation of 4D Semantic City Models

In this section, the second part of the proposed methodology is demonstrated: how to generate 4D semantic datasets from the integrated network of ontologies. To do this, the class and properties assertions (or the TBox) of the ontology network are queried to generate 4D urban data graphs as instances of these assertions (or the ABox). First, these queries must answer the following questions of an atomic datum within a 4D urban dataset:

• Which Class, Property, or Datatype in the network defines the datum? • Does the datum represent geometry?

• Does the datum have temporal properties? SPARQL queries can be implemented to answer these questions. For simplification, a city model as a CityGML 3.0 XML document is considered to be a set of nodes, where each node may have parents and/or children. Determining if a node can be a class or property instance is done using the node's XML tag and searching if a class or property assertion of the same identifier is in the network. Query 1 shows how this can be done for verifying if a node is an object property. This query uses a node's identifier "ntype" and its parent's identifier, "ptype", to find properties with "ptype" in their domain (line 5) or universal quantifications with "ntype" asserted by class restrictions (line 8). A regex function is used to filter these properties with local naming conventions (line 15). Once these types of queries are defined, they can be used to map city model datasets to the ABox of the ontology network as shown in algorithm 1.

Query ] .

14:

} 15:

FILTER regex(STR(?property), "ptype. * .ntype$") 16: } Algorithm 1 looks through all nodes of a given city model in Line 1. For all nodes which are defined by class assertions (line 2) a triple will be created according to if the children of the node (line 4) are defined by a class assertion (line 5), an objectProperty assertion (line 8), a datatypeProperty assertions (line 14), or a datatype assertion (line 21). All of these cases will add the triple (lines 11, 16, 19, 23) to a graph to be returned in line 28. for gc in c.children do During graph generation of CityGML datasets, datatypes from the GeoSPARQL standard are used to represent geometric and geospatial data. To determine when an XML node represents geometry in GML 3.2, GeoSPARQL and stRDF's rules for declaring geo:gmlLiteral values can be implemented according to the definition:

Algorithm 1 Semantic Graph Generation

Valid geo:gmlLiterals are formed by encoding geometry information as a valid element from the GML schema ... in GML 3.2.1 this is every element directly or indirectly indirectly in the substitution group of the element {http://www.opengis.net/ont/gml/3.2}Abstract-Geometry. [START_REF] Ogc | GeoSPARQL -A Geographic Query Language for RDF Data[END_REF] By this definition if an XML node is an instance of a subclass of geo:AbstractGeometry in the ontology network, the XML node represents geometry and its 2D geometric representation can be stored in a triple in the graph. The value of the triple being a copy of the XML data of the original node and its descendant nodes encoded as a literal string. To represent temporal data, xsd:dateTime values are used to represent a point in time. In CityGML 3.0 every city object can use its core:creationDate and core:terminationDate properties as temporal points to representing the beginning and end of its existence.

IMPLEMENTATION AND RESULTING DATASETS

Methodology Implementation

Under the Urban Data Services and Visualization Project (UD-SV) 11 , a transformation tool, UD-Graph 12 , was implemented to 11 https://github.com/VCityTeam/UD-SV 12 https://github.com/VCityTeam/UD-Graph apply this methodology to generate 4D semantic graphs as RDF from CityGML 3.0 data using the RDFLib13 Python library. As shown in figure 4, the methodology proposed in section 3.2 and 3.3 is illustrated. An additional preprocessing step is necessary before graph generation to prepare GML data for use with Geo-SPARQL's 2D queries. This step can add coordinate-referencesystem declarations if necessary and flattens 3D GML geometry to 2D. During generation, gml:id attributes are used as triple identifiers whenever possible. This also applies to Xlinks which are used to link CityGML 3.0's versioning data -such as a city model's version to the corresponding city objects of that version. Afterward, these 4D semantic urban graphs were loaded into a Parliament14 triple store to implement spatio-temporal queries through GeoSPARQL.

The following section will detail how UD-Graph and this methodology are applied to an initial 4D semantic urban dataset.

Resulting Dataset

The 4D semantic urban datasets that are used in these experiments originate from the 1st district of the city of Lyon, France15 as CityGML 2.0. The building and geospatial data from these datasets are converted to CityGML 3.0 with a CityGML 2 to 3 open-source conversion tool 16 and then enriched with temporal versioning data and stored on the CityGML 3.0 Encoding Github17 . Using the proposed transformation methodology, an initial 4D semantic city model dataset was created. Figure 5 illustrates a historical succession of 2 versions of a city model containing 4 buildings. A version transition between these two versions is composed of 3 transactions (or changes): a building deletion (T ransaction 1), a building replacement (T ransaction 2), and a building insertion (T ransaction 3). Figure 6 shows urban data classes and class instances generated from the dataset according to the CityGML building and core ontologies (highlighted blue in figure 1) and the GML geospatial ontologies (highlighted red in figure 1), visualized by OntoGraph 18 .

Query 2 returns all the buildings that have a physical existence intersecting a given point of time. Lines 5-7,9-11 define the creation and termination dates of a building, i.e., the time period of its physical existence. Lines 8 and 12 filter the buildings that exist during the temporal point specified by the user (lines 3,4).

In this example, the query looks for buildings that existed on 01/01/2010. Also, note that the properties in the query use local naming conventions as mentioned in section 3.2.1.

Query 2 Find all buildings that existed at a temporal point ?cityobjectmember a ?type . 7: }

Discussion

This initial approach integrates standardized conceptual models of spatio-temporal and urban data to create integrated 4D semantic datasets as a first step toward 4D semantic city model enrichment. Unlike approaches based on extending the CityGML 2.0 model to support spatio-temporal data, the CityGML 3.0 conceptual model is rich enough to represent 4D semantic city models without extension. By implementing automated conversion tools such as ShapeChange, ISO standard generation of CityGML 3.0 as a ontology network is possible and the benefits of ontology based integration approaches can be taken advantage of with minimal semantic data loss.

The CityGML 3.0 model is especially synergistic with this approach, regarding the creation of 4D semantic data, since the OGC generates the physical data model as XML Schema directly from the UML model using ShapeChange. Therefore, class and properties names from the UML model closely resemble their XML representation and thus transformation mappings from CityGML 3.0 to semantic urban graphs lose minimal semantic data and rest interoperable with their original conceptual model.

In addition, semantic graph formats like RDF align well with CityGML 3.0's graph-like representations of city models in the Versioning module. By integrating geospatial query frameworks like GeoSPARQL spatial navigation of these 4D semantic urban data graphs is possible. This can have many potential applications like change detection of city objects between versions of a city model which may rely on graph formats [START_REF] Nguyen | A Multi-Perspective Approach to Interpreting Spatio-Semantic Changes of Large 3D City Models in CityGML Using a Graph Database[END_REF] and smart city applications based on Semantic Web technologies [START_REF] Gaur | Smart City Architecture and its Applications Based on IoT[END_REF][START_REF] Bischof | Semantic Modelling of Smart City Data[END_REF].

However, there are two modeling limitations of this approach to be addressed in future works. First and foremost, while the CityGML 3.0 application schema is directly generated from the conceptual model, the GML 3.2 model does share this characteristic. Many discrepancies in naming conventions between the GeoSPARQL ontology, ISO 19107 models, and GML 3.2 data were identified. For example, the application schema uses gml:id as a unique local identifier for geometric entities, while this does not exist in the GeoSPARQL ontology and ISO 19107 uses featureID as its identifier. To solve this, additional mappings are being integrated into the ontology network based on ISO TC211 standards to GML 3.2 instance mappings implemented by Enterprise Architect19 , an ISO standard compliant UML modeling platform.

Secondly, as the generated CityGML 3.0 ontologies are constructed from UML's frame-based, "closed-world" assumption of the conceptual model, the mapping transformations to OWL's more "open-world" assumption of the conceptual model are subjective and require some interpretation [START_REF] Brink | Linking spatial data: automated conversion of geo-information models and GML data to RDF[END_REF][START_REF] Cox | An explicit owl representation of iso/ogc observations and measurements[END_REF]. These ontologies are defined according to a more restrictive interpretation in order to guard as many semantic relationships as possible. Because of this, it is possible that while this ontology network works well to generate 4D semantic graphs, other applications of UML to OWL mappings may require a more "open-world" interpretation depending on the purpose of the ontologies.

In addition, to ameliorate the temporal queries proposed in this paper, spatio-temporal query frameworks such as stSPARQL can be implemented to implement more complex temporal queries and improve performance [START_REF] Garbis | Geographica: A benchmark for geospatial rdf stores (long version)[END_REF].

CONCLUSIONS

The evolving urban analysis landscape requires the integration of multiple heterogeneous and autonomous models catering to describing the urban lifecycle. However, a direct translation of data in different evolving data formats to the desired format may lead to semantic data loss. This article presents a methodology for integrating 4D semantic urban data models as a network of ontologies. This UML-based transformation approach based on standards helps to preserve interoperability and reduce semantic loss. This paper also proposes a pipeline for generating integrated RDF data in conformance to this model. The data obtained with this approach were validated and spatio-temporal queries were tested and executed.

Future courses of action include improving the existing semantic model and integrating more urban data models like cadastral data, city documents, and concurrent points of view of urban evolution. Also, the addition of bidirectional transformations could enable RDF to CityGML transformations for non-RDF ready applications. And finally, the implementation of more complex 3D and 4D urban data queries with SPARQL extensions like those provided by stSPARQL and BimSPARQL [START_REF] Zhang | BimSPARQL: Domainspecific functional SPARQL extensions for querying RDF building data[END_REF] are being explored.

DATA AND REPRODUCIBILITY

Detailed notes and the code for reproducing the results of section 3 and 4 can be found on Software Heritage here20 . The generated CityGML 3.0 Versioning dataset is located under the folder /Transformations/test-data/RDF/OWL-basedtransformations.

Figure 1 .

 1 Figure 1. 4D Urban Data Semantic model as a network of ontologies

Figure 2 .

 2 Figure 2. UML Union Flattening example with core:cityModelMember and core:cityObjectMember (top) and their mappings in OWL (bottom)

Figure 3 .

 3 Figure 3. CityGML 3.0 vers:Version UML (top) to OWL Class (bottom)

  Input: A set of ontological entities, O = {o1, o2, ...on}, a set of tree nodes, N = {n1, n2, ...nm} Output: A set of triples { s1, p1, v1 , s2, p2, v2 , ... sn, pn, vn } 1: let graph be an empty set, {} 2: for n in N do 3: if O.hasClass(n) then 4: for c in n.children do 5: if O.hasClass(c) then 6: property = O.findPropertyBetween(n, c) 7: graph.push( n.name, property.name, c.name ) 8: else if O.hasObjectProperty(c) then 9:

Figure 4 .

 4 Figure 4. Proposed 4D Semantic Model Integration and Graph Generation Pipeline

Figure 5 .

 5 Figure 5. Diagram of an RDF Graph of a Historical Succession Transition of Between 2 City Model Versions (above) with several Buildings (below)

  all the city objects in a version, v1, of a concurrent point of view. Line 3 verifies whether the version specified by the user is indeed a version and Line 4 and 5 return all the city object members and type for this given version.Query 3 Find all city objects in a city model version 1: SELECT DISTINCT ?cityobjectmember ?type 2

Figure 6 .

 6 Figure 6. RDF Graph of classes and instances concerning building data generated according to the 4D Urban Data Semantic model defined in Figure 1; CityGML 3.0 Building data is circled in blue; GML geospatial data is circled in red; dashed arrows represent instances of semantic properties (or relationships)

https://www.ogc.org/

https://www.ogc.org/standards/citygml

https://www.isotc211.org/

https://technical.buildingsmart.org/standards/ifc/

https://www.ogc.org/standards/geosparql/

http://www.strabon.di.uoa.gr

https://www.iso.org/standard/57466.html

https://www.w3.org/2001/sw/wiki/Pellet

http://www.hermit-reasoner.com/

https://shapechange.net/

https://rdflib.readthedocs.io/

https://github.com/SemWebCentral/parliament

https://data.grandlyon.com/

https://github.com/tum-gis/citygml2-to-citygml3

https://github.com/opengeospatial/CityGML-3.0Encodings/tree/master/CityGML/Examples

https://sparxsystems.com/enterprise architect user guide/15.2/ model domains/uml profile for gml.html

https://archive.softwareheritage.org/swh:1:dir: e3df50dc7b066ce284ca66c1633e4eb5c2aadf57;origin= https://github.com/VCityTeam/UD-Graph;visit=swh:1:snp: c6ac42c1540a27fbfcdaf09db305d6949751700d;anchor=swh:1: rev:2b02e0c494072dbe4eb6f4d3567c661ae911b2a7

https://projet.liris.cnrs.fr/vcity/

ACKNOWLEDGMENTS

The authors would like to thank University Lumière Lyon 2 for funding the PhD Thesis of Diego Vinasco-Alvarez. This work has been done within the LIRIS Vcity project 21 .