Sarah Benyagoub
email: benyagoub.sarah@univ-mosta.dz

Meriem Ouederni
email: meriem.ouederni@enseeiht.fr

Yamine Ait Ameur
email: yamine@enseeiht.fr

Towards correct Evolution of Conversation Protocols

Keywords: System Evolution, Realisability, Conversation Protocols, Formal Verification, Behavioural Systems

Distributed software systems change dynamically due to the evolution of their environment and/or requirements, their internal designing policies, and/or their specification bugs which must be fixed. Hence, checking system changes must be run continuously. Such systems are usually composed of distributed software entities (called peers) interacting with each other through message exchanges, and this is to fulfil a common goal. The goal is often specified by a conversation protocol (CP), i.e. sequences of sent messages. If there exists a set of peers implementing CP, then CP is said to be realisable. In this paper, we propose a stepwise approach for checking whether an evolution, i.e. adding and/or removing messages and/or peers, can be applied to a CP that was realisable before updating it. We define a set of correct evolution patterns and we suggest an algebra of CP evolution. Our approach ensures that CP evolution preserves the realisability condition.

INTRODUCTION

Distributed software systems change dynamically due to the evolution of their running environment and/or requirements, their internal designing policies, and/or their possible specification bugs which must be fixed. Such systems are usually composed of distributed software entities (called peers) evolving concurrently in a distributed setting and interacting with each other throughout messages exchanges to fulfil a common goal.

In a top-down design of distributed software, the interaction among peers is usually modelled using collaboration diagrams, Message Sequence Chains (MSCs) or conversation protocols (CP) [START_REF] Bultan | Modeling interactions of web software[END_REF]). Let us focus on CPs, these describe interactions among peers by describing uniquely the allowed sequences of sent messages. Here, it is crucial to know if the set of interactions in a CP can be implemented. In other word, considering a CP, one must check whether there exists a set of peers where their composition generates the same sequences of send messages as specified by the CP. This issue characterises the realisability problem.

In order to formally specify, verify, and fix issues which violate realisability, CP can be modelled using Labelled Transition Systems (LTSs) where communication follows either synchronous or asynchronous semantics. Using this model enables automated analysis of interaction properties, e.g., realisability checking. Although it is obvious to check realisability in the case of synchronous communication, this realisability problem remains undecidable in the most general setting [START_REF] Brand | On Communicating Finite-State Machines[END_REF] due the possibly ever-increasing queuing mechanism and unbounded buffers in the case of asynchronous communication.

The recent work of [START_REF] Basu | Deciding Choreography Realizability[END_REF]) proposed a necessary and sufficient condition for verifying if a CP can be implemented by a set of peers communicating asynchronously throughout FIFO buffers with no restriction on their buffer sizes. This work solves the realisability issue for a subclass of asynchronously communicating peers, namely, the synchronisable systems, i.e., the system composed of interacting peers behaves equally by applying synchronous or asynchronous communication. A CP is realisable if there exists a set of peers implementing that CP, i.e., they send messages to each other in the same order as CP, and such that their composition is synchronisable. In [START_REF] Basu | Deciding Choreography Realizability[END_REF], the full checking of CP realisability applies the following steps: i) peer projection from CP; ii) checking synchronisability; and iii) checking equivalence between CP and its distributed system.

Based on LTS model, we can verify CP realisability using existing techniques such as model checking for systems with reasonable sizes [START_REF] Basu | Deciding Choreography Realizability[END_REF]) (i.e., number of states, transitions and communicating peers) or theorem proving for scalable systems [START_REF] Farah | A Correct-by-Construction Model for Asynchronously Communicating Systems[END_REF].

Considering realisable CPs, we are interested in studying the evolution of those CPs. In fact, these specify cross-organisational interactions with no centralised control between peers which can be administrated and executed by geographically distributed and autonomous companies. In order to cope with new environment changes and enduser requirements, system interaction and the corresponding CP need to evolve continuously over time. However, changing CP might result in knockon effects on its realisability. Hence, verifying the correctness of CP evolution to ensure realisability preservation must also be run continuously.

Regarding the literature, existing work such as (Rinderle et al. 2006b;[START_REF] Ryu | Supporting the Dynamic Evolution of Web Service Protocols in Serviceoriented Architectures[END_REF][START_REF] Roohi | Realizability and Dynamic Reconfiguration of Chor Specifications[END_REF] give some solutions for system evolution. In (Rinderle et al. 2006b;[START_REF] Ryu | Supporting the Dynamic Evolution of Web Service Protocols in Serviceoriented Architectures[END_REF] the authors propagate the choreography updates into communicating peers. [START_REF] Roohi | Realizability and Dynamic Reconfiguration of Chor Specifications[END_REF] focuses on system reconfiguration meaning that for a CP which has been updated into CP', the authors check whether the traces that has been executed in CP can be performed again in CP'. This reconfiguration can be better applied for run-time system to ensure execution consistency. All these approaches do not consider realisability preservation.

There exist other research approaches which can be applied as a posteriori evolution checking. The approaches suggest solutions every time the realisability check fails. For example, existing work on enforcing CP realisability such as the one given in (G üdemann et al. 2012) and recently on CP repairability [START_REF] Basu | Automated Choreography Repair[END_REF]) can be used to ensure the realisability of an already updated CP.

Our statement is different than existing work and it is as follows: an evolution is allowed if it does not violate the CP realisability. By doing so, we suggest a priori verification approach of CP evolution. Instead of running the full realisability checking as described previously and detailed in Section 2, our proposal consists in performing partial verification uniquely at the CP level in order to answer the question if there still exist a set of peers implementing the updated CP. In this work, we consider the evolution at the CP level and we study its realisability effect on the distributed peers. The main issue is considering that system specifications may change over time (e.g., service upgrade or degrade by adding and/or removing either messages exchanges or interacting peers), how can we ensure realisability preservation? to answer these questions, we proceed as follows:

• We first describe CP using LTS

• We rely on the realisability condition given in [START_REF] Basu | Deciding Choreography Realizability[END_REF] • We identify the set of behavioural properties which can hold by CP evolution yet they violate the realisability condition

• We suggest a set of evolution patterns and we show how the application of such patterns do not violate CP realisability

• We propose a language for correct CP evolution

The remainder of this paper is structured as follows: Section 2 introduces the formal definitions and the background on which our proposal relies. Section 3 presents the behavioural properties to be checked before application of CP evolution. In Section 4, we suggest an algebra for CP evolution with no violation of realisability. We present in Section 5 a case study to illustrate our approach. Section 6 overviews related work. Lastly, we conclude our work and present some perspectives in Section 7.

DEFINITIONS

In this section, we present our behavioural model for peers and CP. We, then, define how distributed peers can be obtained by projection from a given CP. Lastly, we define synchronisable systems, and we present realisability condition considering asynchronous communication.

Peer Model

We use Labeled Transition Systems (LTSs) for modelling the CP and the peers included in that specification. This behavioural model defines the order of sent messages in CP. At the peers level, the LTS can be computed by projection from CP and they describe the order in which those peers execute their send and receive actions.

Definition 1 (Peer) A peer is an LTS P = (S, s 0 , Σ, T ) where S is a finite set of states, s 0 ∈ S is the initial state, Σ = Σ ! ∪ Σ ? ∪ {τ } is a finite alphabet partitioned into a set of send messages, receive messages, and the internal action, and T ⊆ S ×Σ×S is a transition relation.

We write m! for a send message m ∈ Σ ! and m? for a receive message m ∈ Σ ? . We use the symbol τ (tau in figures) for representing internal activities. A transition is represented as s l → s where l ∈ Σ.

(internal) s τ → s ∈ T a if ∃i ∈ {1, . . . , n}, (i) s i τ → s i ∈ T i , (ii) ∀k ∈ {1, . . . , n} : Q k = Q k , and (iii) ∀k ∈ {1, . . . , n} : k = i ⇒ s k = s k
We use LT S k a to define the bounded asynchronous composition, where each message buffer is bounded to size k. The definition of LT S k a can be obtained from Def. 5 such that the maximum number of send messages that can be stored in the buffers is limited to k. A system is synchronizable [START_REF] Basu | Deciding Choreography Realizability[END_REF] when its behavior remains the same under both synchronous and asynchronous communication semantics. This is checked by bounding buffers to k = 1 and comparing interactions in the synchronous system with send messages in the asynchronous system.

Definition 6 (Synchronizability) Given a set of peers {P 1 , . . . , P n }, the synchronous system (P 1 | . . . | P n ) = (S s , s 0 s , L s , T s ), and the 1-bounded asynchronous system

((P 1 , Q 1 1 ) || . . . || (P n , Q 1 n )) = (S a , s 0
a , L a , T a ), two states r ∈ S s and s ∈ S a are synchronizable if there exists a relation Sync such that Sync(r, s) and:

• for each r m -→ r ∈ T s , there exists s m! --→ s ∈ T a , such that Sync(r , s ); • for each s m! --→ s ∈ T a , there exists r m -→ r ∈ T s , such that Sync(r , s ); • for each s m? --→ s ∈ T a , Sync(r, s ).
The set of peers is synchronizable if Sync(s 0 s , s 0 a ).

Example 4

The system described in Figure 1 is not synchronisable because peer 1 can send "a" and "c" in sequence in the asynchronous system, whereas "b" occurs before "c" in the synchronous system as specified in the CP.

In order to check CP realisability, there is a need to check well-formedness. Well-formedness states that whenever the i-th peer buffer Q i is nonempty, the system can eventually move to a state where Q i is empty. For every synchronizable set of peers, if the peers are deterministic, i.e., for every state, the possible send messages are unique, wellformedness is implied.

The approach presented in [START_REF] Basu | Deciding Choreography Realizability[END_REF] proposes a sufficient and necessary condition showing that the realizability of conversation protocols is decidable.

Definition 7 (Realizability) A conversation protocol CP is realizable if and only if (i) the peers computed by projection from this protocol are synchronizable, (ii) the 1-bounded system resulting from the peer composition is well-formed, and (iii) the synchronous version of the distributed system {P 1 , . . . , P n } is equivalent to CP .

In the remainder of this paper, we refer to a realisable CP as R(CP).

Both synchronizability and realizability properties can be checked automatically using equivalence checking as done in [START_REF] Basu | Deciding Choreography Realizability[END_REF]. This check requires the modification of the asynchronous system for hiding receptions (m? τ ), renaming emissions into interactions (m! m), and removing τtransitions using standard minimization techniques. The correctness of that approach is given in [START_REF] Farah | A Correct-by-Construction Model for Asynchronously Communicating Systems[END_REF].

BEHAVIOURAL PROPERTIES

In order to check the realisability of a CP that has been updated, we must ensure that the resulting LTS does not hold some branching and/or sequential structures which violate realisability condition. We define in the following some properties which enable us to check such structures.

Branches related Properties Property 1 (Non-Deterministic Choice (NDC))

Given a conversation protocol CP =< S CP , s 0 CP , L CP , T CP >, a state s CP ∈ S CP is called non-deterministic branching state if :

• ∃{s CP m P i ,P j -----→ s CP , s CP m P i ,P j -----→ s CP } ⊆ T CP , and

• s CP = s CP
This choice is referred to as non-deterministic choice.

We define in the following divergent choice (this is also called non-local branching choice in the literature) and it is different than process divergence (Ben-Abdallah and Leue 1997).

Property 2 (Divergent-Choice) Given a conversation protocol CP =< S CP , s 0 CP , L CP , T CP >, a state s CP ∈ S CP is divergent branching state if :

• ∃{s

m P i ,P j -----→ s CP , s CP m P j ,P i -----→ s CP } ⊆ T CP ,

and

• s CP = s CP , and

• m = m
This choice is referred to as divergent choice.

Sequences related Properties

Given a CP, there is at least two partitions of peers where there exist no interaction between both partitions.

Property 3 (Independent Sequences (ISeq))

Given a conversation protocol CP =< S CP , s 0 CP , L CP , T CP >, a transition sequence

s CP m P i ,P j -----→ . . . s CP m P k ,Pq
-----→ s CP , where all transitions are in T CP , is called independent sequence if:

• {P i , P j } ∩ {P k , P q } = ∅
The following property enables us to detect sequences in CP which lead to non-local emission choices made by two different peers in the distributed system. To avoid that situation, every peer that join the conversation at an intermediate state (i.e., different than the initial state) must be receiver the first time it shows up. Otherwise, if a peer would like to send a message m at an intermediate state, then this must be receiver in its last interaction before sending m.

Property 4 (Divergent Sequences (DSeq))

Given a CP, there exists a transition sequence

s 0 CP m P i ,P j -----→ . . . s CP m P k ,Pq
-----→ s CP where all transitions are in T CP :

• for every sender peer P t appearing before state s CP , t = k, or

• there is at least a transition s CP m P k ,P t -----→ s CP ∈ T CP such that:

s CP is reachable from s CP , and

-there is no transition in s CP m P k ,P t -----→ s CP . . . s m P k ,Pq
-----→ s CP where P k is receiver.

COMPOSITIONAL REALISABILITY

CP evolution stands for two possible tasks, namely, adding and/or removing either messages and/or interacting peers. We define here how CP realisability can be preserved by applying some evolution patterns presented in the following.

Evolution Patterns

We introduce in this paper two composition operators denoted as ⊗ (+,sCP ) for branching composition and ⊗ ( ,sCP ) for sequential composition. We also assume other operators not presented here for lack of space, namely, ⊗ ( ,sCP ) for parallel composition, and ⊗ ( ,sCP ) for looping composition. The operator ⊗ ( ,sCP ) generates at a state s CP all the interleaved behaviour of a set of transitions such that every generated branch must satisfy sequence related properties. The operator ⊗ ( ,sCP ) enables us to add self-loop of the form s m P i ,P j -----→ s where i = j and such that sequence related properties must be preserved. , L CP i , T CP i > and a state s CP ∈ S CP , the branching composition ⊗ (+,sCP ) (CP, {CP 1 , . . . , CP n }) means that there is a choice at s CP between the remaining behaviour of CP (i.e., starting from s CP ) and all CP i such that:

• Properties 1 and 2 do not hold at the state s CP , and

• ∀CP i , ⊗ ( ,sCP ) (CP, CP i ) holds
Remark 1 The application of an operator ⊗ (op,sCP ) (CP , CP ) assumes that the initial state of CP is fused with the state s CP .

Algebra for CP Evolution: syntax and language

We introduce in Listing 1 a CP algebra which we use for defining the evolution such that realisability is preserved. We refer to a state s f as final if there is no outgoing transition at that state. We denote by ECP a CP that evolutes over time while preserving realisability. The expression ECP + stands for one or more ECP .

ECP ::= ECP b | ECP op ECP b + ECP b ::= s (Pi,m,Pj ) ------→ s | ∅ op ::= ⊗ (+,s f ) | ⊗ ( ,s f ) | ⊗ ( ,s f ) | ⊗ ( ,s f )
Listing 1: CP Evolution Grammar

Realisable CP Evolution

Conjecture 1 ECP b is realisable.

Proof 1 This is obvious by default.

Conjecture 2 Given an ECP =< S ECP , s 0 ECP , L ECP , T ECP > and a ECP b such that R(ECP ) and R(ECP b ), s f ∈ S ECP , then requires to check its effect on other partner peers. This approach is implemented into DYCHOR.

In [START_REF] Preda | Dynamic Choreographies -Safe Runtime Updates of Distributed Applications[END_REF][START_REF] Fdhila | Dealing with change in process choreographies: Design and implementation of propagation algorithms[END_REF], both approaches study the evolution that might arise at the peers side. The authors propagate the change from one peer to other partners. The work proposed in [START_REF] Fdhila | Dealing with change in process choreographies: Design and implementation of propagation algorithms[END_REF] applies to Business Process Management (BPM) domain and Service Oriented Architecture (SOA). It describes service choreographies using tree-based model. The authors consider some changes such as "replace, delete, update, and insert" of behavioural fragments. [START_REF] Preda | Dynamic Choreographies -Safe Runtime Updates of Distributed Applications[END_REF] define a new language referred to as DIOC for programming distributed applications that are free from deadlocks and races by construction. The semantics of DIOC language relies on labelled transition systems. The approach given in [START_REF] Preda | Dynamic Choreographies -Safe Runtime Updates of Distributed Applications[END_REF] enables to update fragment of codes of distributed peers. This can be specified at the choreography level where blocks of code that can be updated dynamically must be tagged using "scope". These "scope" blocks have to be known a priori when describing the choreography. The solutions given in [START_REF] Preda | Dynamic Choreographies -Safe Runtime Updates of Distributed Applications[END_REF][START_REF] Roohi | Realizability and Dynamic Reconfiguration of Chor Specifications[END_REF][START_REF] Jureta | Dynamic Requirements Specification for adaptable and open Service-oriented Systems[END_REF] deal with run-time evolution.

To the best of our knowledge, we are the first who verify the evolution at the CP level such that its realisability must be preserved. Furthermore, we have no restriction on the application domain, yet we use a formal model which can be applied for design, verification, and implementation of different distributed systems, e.g. Web services, concurrent systems, Cyber Physical Systems, etc. Our result applies also for asynchronously communicating systems as far as these are synchronisable with no restriction on the buffer size.

CONCLUSION AND PERSPECTIVES

In this paper, we presented a preliminary solution for correct evolution of distributed system for which their interaction is described with a conversation protocol. We proposed a language which enables one to incrementally design distributed system that can be updated over time such that their realisability is preserved while applying and composing evolution operators.

In the near future we aim at formalising all properties and composition operators used in this paper. We intend to prove that our evolution operators and their properties preserve CP realisability. Our conjectures have also to be formally proved. We aim also at defining looping and parallel operators as well as extending our language with new operators for messages broadcast and multicast. Lastly, we are using theorem proving techniques in order to prove the correctness of CP evolution. Based on proofbased techniques, we aim at handling any number of peers and exchanged messages such that scalability is ensured.

  CP , T CP > and a state s CP ∈ S CP , the sequential composition ⊗ ( ,sCP ) (CP, CP ) means that CP must be executed before CP such that Properties 3 and 4 do not hold. CP , T CP >, a set {CP i }, i = 1..n such that CP i =< S CP i , s 0

	Definition 8 ⊗ ( ,sCP )	Given		a	CP	=<
	S CP , s 0 CP , L CP , T CP S CP , s 0 CP , L Definition 9 ⊗ (+,sCP )	>, Given	a	a	CP CP	=< =<
	S CP , s 0 CP , L CP i


Théo Berthet

Jean-François Ritz

Matthieu Ferry

Phuntsho Pelgay

Rodolphe Cattin

Dowchu Drukpa

Regis Braucher

György Hetényi

Active tectonics of the eastern Himalaya: New constraints from the first tectonic geomorphology study in southern Bhutan

de 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 AND TECTONIC SETTING

The Himalayan belt is one of the few on-land collision zones where great earthquakes with sizes comparable to those of oceanic subduction zones can occur. Over the past few centuries several major earthquakes (M > 8) have been documented both geologically and historically in the central and eastern Himalaya: in western Nepal (A.D. 1505, ~M8.7), in Bihar-Nepal (1934, M8.1), in Assam (1950, ~M8.4), and near the Shillong plateau (1897, ~M8.1) [START_REF] Kumar | Paleoseismological evidence of surface faulting along the northeastern Himalayan front, India: Timing, size, and spatial extent of great earthquakes[END_REF][START_REF] Lavé | Evidence for a great medieval earthquake (1100 A.D.) in the central Himalayas[END_REF][START_REF] Yule | Possible evidence for surface rupture of the Main Frontal Thrust during the great 1505 Himalayan earthquake, far-western Nepal[END_REF][START_REF] Mugnier | Structural interpretation of the great earthquakes of the last millennium in the central Himalaya[END_REF] (Fig. 1). In this context, the Bhutan Himalaya, which is located between the great Himalayan ruptures of 1934 and 1950, appears as an ~350-km-long section where no similar great earthquakes have ever been documented. A single historical account reports an earthquake in 1713, but its magnitude and accurate location are unknown [START_REF] Ambraseys | A note on early earthquakes in northern India and southern Tibet[END_REF]. The occurrence of the 1897 event in the Shillong plateau attests to the accommodation of some amount of convergence in this region, and has been proposed to increase the interval between great earthquakes in the Bhutan Himalaya [START_REF] Bilham | Plateau 'pop-up' in the great 1897 Assam earthquake[END_REF]. [START_REF] Gahalaut | Low seismicity in the Bhutan Himalaya and the stress shadow of the 1897 Shillong Plateau earthquake[END_REF] also proposed that the stress shadow caused by this earthquake may be responsible for the low seismicity rate currently observed in Bhutan.

In terms of paleoseismicity, a study by [START_REF] Lavé | Evidence for a great medieval earthquake (1100 A.D.) in the central Himalayas[END_REF] along the Main Frontal thrust (MFT) in central Nepal suggests the occurrence of a great medieval earthquake at ca. A.D. 1100, with an estimated vertical slip component of 7-7.5 m. From similar investigations east and west of Bhutan, [START_REF] Kumar | Paleoseismological evidence of surface faulting along the northeastern Himalayan front, India: Timing, size, and spatial extent of great earthquakes[END_REF] proposed that this medieval event may have broken an 800-km-long portion of the MFT, including Bhutan, with a magnitude approaching M9. However, a more recent study by [START_REF] Sapkota | Primary surface ruptures of the great Himalayan earthquakes in 1934 and 1255[END_REF] in eastern Nepal documented two great earthquake ruptures contemporary with the 1934 and 1255 historical earthquakes, instead of a single, giant 11 th century event.

In central Nepal and Arunachal Pradesh (east of Bhutan), tectonic geomorphology studies have determined Holocene slip rates along the MFT of 21 ± 1.5 mm/yr [START_REF] Lavé | Active folding of fl uvial terraces across the Siwaliks Hills, Himalayas of central Nepal[END_REF] and 23 ± 6.2 mm/yr [START_REF] Burgess | Holocene shortening across the Main Frontal Thrust zone in the eastern Himalaya[END_REF], respectively. In Nepal, a consistent shortening rate of ~18-20 mm/yr was measured by Active tectonics of the eastern Himalaya: New constraints from the fi rst tectonic geomorphology study in southern Bhutan GPS (e.g., [START_REF] Ader | Convergence rate across the Nepal Himalaya and interseismic coupling on the Main Himalayan Thrust: Implications for seismic hazard[END_REF] (Fig. 1). This suggests that the interseismic loading accumulated beneath the High Himalaya on the Main Himalayan thrust (MHT, the down-dip continuation of the MFT and the main interface along which India underthrusts the Himalaya) is relaxed during large events that transfer the shortening to the emergent MFT [START_REF] Cattin | Modeling mountain building and the seismic cycle in the Himalaya of Nepal[END_REF][START_REF] Mugnier | Structural interpretation of the great earthquakes of the last millennium in the central Himalaya[END_REF].

The peculiar tectonic setting of the Bhutan Himalaya raises the question of whether its seismic behavior differs from what is occurring in Nepal. To answer this question, we carried out a tectonic geomorphology study in south-central Bhutan that yields the fi rst estimate of the Holocene slip rate along the Bhutanese frontal thrust, and the fi rst constraints on major seismic events that occurred during the past millennium.

TECTONIC GEOMORPHOLOGY ANALYSIS

We focused our work in southern Bhutan at 90°-90.5°E within the Sarpang reentrant (Fig. 2A) where the topographic front is 15 km farther north compared to the rest of the region (Long et al., 2011). South of the topographic front, a ridge of hills rises above the surrounding alluvial fans. It spans 30 km laterally and ends at longitude 90.3°E (Fig. 2A). This structure uplifts and defl ects the drainage network and illustrates the propagation of the deformation to the south [START_REF] Yeats | Active faulting south of the Himalayan Front: Establishing a new plate boundary[END_REF]. From the analysis of a 5-m-resolution digital elevation model (DEM) extracted from Pleiades satellite images and 1:25,000-scale aerial stereoscopic photographs, we selected two sites along the topographic front at the outlets of the Sarpang Khola and the Mao Khola rivers (Fig. 2B). In both areas, tectonic scarps and uplifted alluvial terraces attest to the accumulation of vertical deformation through time. Note that herein, terraces are named incrementally starting with the surface topographically closest to the active stream. The fi rst site is located on the western part of the reentrant within the town of Sarpang; the second one on the eastern part north of the town of Gelephu. In the Sarpang area (Fig. 2B), abandoned terraces are preserved on both banks of the river. Two main terraces (T1 and T2) can be distinguished, and are affected to the south by an east-west-trending fault scarp perpendicular to the drainage. In the Gelephu area (Fig. 2C), we observe similar fl uvial terraces, T1 and T2, plus an additional higher one named T3 preserved on the east bank of the river. Within the west bank, the east-west-trending scarp bounds the uplifted terraces T1 and T2. Within the east bank of the river, the scarp is found 1.5 km farther north and trends also east-west, cutting through the T1 terrace. There we assume that the two scarp segments are connected through a northeast-southwest section within the river bed. Further east, the trend of the scarp affecting T1 veers abruptly north-south and splits into two fault strands within T2 and T3 on the east bank of the river. Along that section, deformation is partitioned between a reverse fault located at the foot of the uplifted terraces and a right-lateral strike-slip fault forming the boundary between T2 and T3 (a similar feature was observed at the Kala Amb site of [START_REF] Kumar | Paleoseismological evidence of surface faulting along the northeastern Himalayan front, India: Timing, size, and spatial extent of great earthquakes[END_REF]).

Sarpang Area

Figure 3A shows a perspective view of the east bank of the Sarpang Khola where we focus our study. The lower and youngest affected terrace, T1, is preserved along the riverbed of the Sarpang Khola (Fig. 3A). Figure 3B shows the riser of terrace T1 made of cobble and boulder layers topped with a sand unit. The terrace is affected by an ~4-m-high fault scarp perpendicular to the fl ow of the river. Within the hanging wall, immediately north of the scarp, one observes the slightly folded strath horizon of terrace T1 above a highly fractured bedrock unit (Baxa Forma-tion). These features suggest that the Baxa Formation overthrusts the Quaternary deposits.

Terrace T2 is delimited to the south by an abrupt east-west topographic front perpendicular to the fl ow of the Sarpang Khola, defi ning a cumulative fault scarp of several tens of meters in height that corresponds to the eastward extension of the east-west-trending scarp affecting terrace T1. Topographic profi les of T2 across this scarp from a kinematic GPS survey and the 5-m-resolution DEM allow estimating a cumulative offset of 53.4 ± 2.2 m (Fig. 3C). This corresponds to a minimum tectonic offset, as the footwall can be partially buried by colluvium and alluvium; this applies to all tectonic offsets measured in this study. At the eastern foot of the cumulative scarp we observe a smaller scarp of 4.4 ± 0.5 m, comparable to the one affecting terrace T1 (Fig. 3D).

The abandonment age of the fl uvial terrace T2 was determined by in situ-produced cosmogenic 10 Be from the measurement of 10 Be concentrations along a depth profi le. Three stratigraphic units are observed in the pit: (1) an upper 0.5-mthick sandy unit, (2) an intermediate clayey sand unit with angular pebbles (colluvium) radiocarbon-dated at A.D. 795 ± 105, and (3) a lower 1.5-m-thick sand unit with decimetric cobbles (debris fl ow). The colluvial nature of the intermediate layer and its young age attest to the burial of the underlying terrace. To model 10 Be concentrations at depth, we account for burial since A.D. 795 ± 105 and obtain an exposure age of 6400 ± 1300 yr (2σ uncertainty) for terrace T2 (see details in the GSA Data Repository1 ).

To date the deposition of terrace T1, we collected charcoal samples within a 70-cm-deep pit dug into the sand unit found on top of terrace T1 within the hanging wall. Two sub-units can be distinguished. Two charcoals within the upper sub-unit yield modern dates and are interpreted as being anthropic (the terrace surface has been reworked for construction purposes). The two charcoals collected within the lower sub-unit give consistent dates of A.D. 1520 ± 100 and A.D. 1570 ± 80. We interpret them as corresponding to the depositional age of terrace T1 (see details in the Data Repository).

Gelephu Area

Figure 4 focuses on our study area on the east bank of the Mao Khola, near Gelephu. At this site, the terrace T1 can be subdivided into two sub-terraces (T1a and T1b) that are both affected by an ENE-WSW-trending fault scarp. Note that T1a and T1b are found on both sides of the scarp (due to river incision). We surveyed this area using kinematic GPS to build a highresolution DEM with high vertical accuracy (Fig. 4A). The map shows that the fault scarp is buried under a debris fan where its strike changes from ENE-WSW to NW-SE. The vertical offset of terrace T1a is 4.6 ± 0.5 m (Fig. 4B). The vertical offset of terrace T1b is 7.9 ± 0.5 m (Fig. 4C) across the ENE-WSW scarp, and 7.4 ± 0.5 m (Fig. 4D) across the NW-SE scarp.

We used radiocarbon and in situ-produced 10

Be to date the deposition of terrace T1b. The

Holocene Slip Rate along the Bhutanese Frontal Thrust

Combining the mean offset of 53.4 ± 2.2 m across terrace T2 in Sarpang and its 10 Be age of 6400 ± 1300 yr yields a mean vertical slip rate of 8.8 ± 2.1 mm/yr along the Bhutanese frontal thrust for the Holocene period (Fig. 5). Offsets and ages obtained near Sarpang on terrace T1 and near Gelephu on terrace T1b are consistent with this estimate (Fig. 5). From this vertical slip rate, we estimate the slip rate along the thrust fault: we calculate a dip angle of ~20° for the thrust at the surface using the horizontal shift of the fault trace between the terraces T1a and T1b at the Gelephu site. However, this value may be an underestimate because of lateral incision at the foot of the scarp affecting T1a. At the Sarpang site, where we calculated the vertical slip rate, Long et al. (2011) reported structural measurements of dip bedding of ~30°. Considering both values, we use a dip of 25° ± 5° together with mean vertical slip rate of 8.8 ± 2.1 mm/yr to calculate a Holocene slip rate along the fault of 20.8 ± 8.8 mm/yr (see details in the Data Repository).

Age Constraints on the Two Last Great Earthquakes in Bhutan

Our observations in the Sarpang and Gelephu areas show evidence for an ~4.5-m-high fault scarp affecting the most recent alluvial terrace T1 (see Figs. 3B, 3D, and4B). We interpret this offset as the result of the latest surface-rupturing event along the Bhutanese frontal thrust. If we consider a 20° dip angle for this thrust at the surface, this vertical displacement would result from an ~13 m slip along dip. The radiocarbon age obtained for the deposition of terrace T1 in Sarpang indicates that the fi rst event occurred after A.D. 1570 ± 80. This event may correspond to the A.D. 1713 event, the only historical earthquake reported for the region [START_REF] Ambraseys | A note on early earthquakes in northern India and southern Tibet[END_REF]. At Gelephu, T1b is affected by a fault scarp of ~8 m, which is about twice the height of the scarp affecting T1a (4.6 m). We consider that T1b has been displaced by two events with vertical components of ~4.5 m for the ultimate and of ~3.5 m for the penultimate. This yields an ~10 m slip for the penultimate event considering a 20° dip for the fault. Age constraints on T1b at Gelephu indicate that the penultimate event occurred after A.D. 1150 ± 100.

The displacements observed for these two Bhutanese events are similar to co-seismic displacements associated with the M8.1 1934 earthquake in eastern Nepal [START_REF] Sapkota | Primary surface ruptures of the great Himalayan earthquakes in 1934 and 1255[END_REF], suggesting the two last surface-rupturing events along the Bhutanese frontal thrust were great earthquakes with magnitudes that may have been above M8.

DISCUSSION AND CONCLUSIONS

Our study across the Bhutanese frontal thrust brings new constraints on the active tectonics in the eastern Himalaya. A Holocene vertical slip rate of 8.8 ± 2.1 mm/yr was estimated. This yields a horizontal shortening rate of 20.8 ± 8.8 mm/yr considering a dip of 25° ± 5°. This estimate is consistent with shortening rates of 21.5 ± 1.5 mm/yr and 23.0 ± 6.2 mm/yr estimated in central Nepal [START_REF] Lavé | Active folding of fl uvial terraces across the Siwaliks Hills, Himalayas of central Nepal[END_REF] and in India east of Bhutan [START_REF] Burgess | Holocene shortening across the Main Frontal Thrust zone in the eastern Himalaya[END_REF], respectively. This suggests that the Himalayan convergence is also mainly accommodated on the MFT at the longitude of Bhutan. [START_REF] Drukpa | GPS constraints on Indo-Asian convergence in the Bhutan Himalaya: Segmentation and potential for a 8.2<Mw<8.8 earthquake[END_REF] gave a preliminary slip rate estimate of ~20 mm/yr from new GPS observations, consistent with the slip rate estimated in this study. This estimate suggests that the deformation accommodated within the foreberg located in the Indian plain south of Sarpang might be small. Furthermore, our observations document two major events with M > 8, cumulating ~8 m of vertical deformation, that have occurred since A.D. 1150 ± 100. The most recent event occurred after A.D. 1570 ± 80 and could correspond to the A.D. 1713 historical event. If so, the poorly constrained magnitude (Mw ~7) proposed by [START_REF] Ambraseys | A note on early earthquakes in northern India and southern Tibet[END_REF] for the 1713 event may be an underestimate. This is not inconsistent with Jackson's (2002, p. 147) report mentioning that "occurred [in Bhutan] the great terror of an earthquake that pulverized all houses and huts in every direction." The penultimate event occurred between A.D. 1150 ± 100 and A.D. 1570 ± 80. This age constraint is not inconsistent with a potential giant medieval earthquake at ca. A.D. 1100. However, the occurrence of two surfacerupturing events during the past millennium in Bhutan is similar to observations made by [START_REF] Sapkota | Primary surface ruptures of the great Himalayan earthquakes in 1934 and 1255[END_REF] in eastern Nepal where two major earthquakes occurred in A.D. 1934 and A.D. 1255. We thus favor the interpretation of strain being accommodated into individual rupture events, along with the role of fault segmentation in limiting rupture size.

To conclude, the fact that a similar amount of Holocene deformation is accommodated along the MFT in Bhutan as in the neighboring parts of the arc suggests that there is no major lateral variation of the distribution of strain along the frontal thrust from central Nepal to east of Bhutan despite the presence of the Shillong plateau. Moreover, we confi rm that most of the strain accumulated along this major convergent system is accumulated along the frontal thrust. The two surface ruptures observed in this study also show that the Bhutan Himalaya is able to generate great Himalayan events and that the continental shortening processes beneath the Shillong plateau have little effect in reducing their propensity. Thus, the size of rupture events in the Himalaya seems most likely related to the role of fault segmentation rather than major longitudinal variations of convergence rate.

3

  Aix-Marseille Université, CNRS-IRD-Collège de France, CEREGE UM34, Aix-en-Provence 13545, France 4 Swiss Seismological Service, ETH-Zürich, 8092 Zurich, Switzerland 5 Department of Earth Sciences, ETH-Zürich, 8092 Zurich, Switzerland ABSTRACT How convergent systems distribute strain among frontal thrusts is a major concern regarding seismic hazard assessment. Along the 2500 km Himalayan arc, the seismic behavior of the Bhutan region is unknown, because it corresponds to the only portion of the arc where no evidence of major earthquakes has been reported. This can be due either to the fact that no active tectonic studies have been conducted or to continental shortening being absorbed by the Shillong plateau 150 km farther south. Analyzing offset fl uvial terraces in south-central Bhutan shows that two major earthquakes ruptured the Himalayan frontal thrust during the last millennium, and that a comparable rate of Holocene deformation (~20 mm/yr) is accommodated across the Himalaya in Bhutan as in central Nepal. Thus, the propensity for great earthquakes in Bhutan is similar to what is observed in neighboring portions of the Himalaya arc. This in turn suggests that the shortening process beneath the Shillong plateau has little effect on how strain accumulates within the Bhutanese Himalaya.

  Fig 2A

Figure

  Figure 2. A: Shaded SRTM3 (Shuttle Radar Topography Mission) topography of Sarpang reentrant in central Bhutan with drainage (blue) and active thrust faults (red). B, C: Detailed topographic maps of Sarpang and Gelephu sites showing uplifted river terraces. Dashed red line corresponds to our interpreted fault trace.

2

  -m-deep soil pit dug into the terrace surface shows two main depositional units, both clast supported and topped by sandy-silty horizons. A brown organic-poor modern soil has developed in the silty upper part of the fi rst (upper) unit. Two radiocarbon dates obtained from a bulk sample collected in the sands that topped the second (lower) unit yield close dates of A.D. 955 ± 65 and A.D. 1150 ± 100. Considering the oldest dated sample (A.D. 955 ± 65) as reworked, we retain A.D. 1150 ± 100 as the date constraint for T1b (see details in the Data Repository).

  Figure 3. A: Perspective view to northeast of Sarpang Khola east bank (Pleiades satellite image). Red and orange lines defi ne trace of fault and profi les across terrace T2, respectively. B: Fault scarp affecting T1 terrace. C: Topographic profi les across T2. D: Photo and topographic profi le of fault scarp at foot of T2.

Figure 4 .

 4 Figure 4. A: Map of T1 terrace fl ight observed on Mao Khola east bank (see location in Fig. 2C). Contour lines are every 2 m; red line corresponds to fault trace. B, C, D: Profi les and photos of fault scarps affecting T1a and T1b (see black lines on A for their locations). Note debris fan burying scarp.

Figure

  Figure 5. Vertical slip rate along Bhutanese frontal thrust (black line) and its uncertainties (dashed lines). Rectangles represent error on ages and elevations of terraces, respectively. Arrows within T1a and T1b rectangles indicate that the seismic event affecting each terrace is younger than their respective depositional age.

  Théo Berthet 1 , Jean-François Ritz 1 , Matthieu Ferry 1 , Phuntsho Pelgay 2 , Rodolphe Cattin 1 , Dowchu Drukpa 2 , Régis Braucher 3 , and György Hetényi 4,5 1 Géosciences Montpellier, UMR5243, Université Montpellier 2, Place E. Bataillon, 34095 Montpellier, France 2 Seismology and Geophysics Division, Department of Geology and Mines, Ministry of Economic Affairs, Post Box 173, Thimphu, Bhutan

GSA Data Repository item

2014152, details on dating and uncertainties, is available online at www .geosociety.org/pubs/ft2014.htm, or on request from editing@geosociety.org or Documents Secretary, GSA, P.O. Box 9140, Boulder, CO 80301, USA.

ACKNOWLEDGMENTS

We are grateful to R. Bilham, J. Hollingsworth, and two anonymous reviewers for their helpful comments and suggestions. We thank K. Thinley (Land Survey of Bhutan) for providing aerial photographs and the AS-TER Team for their assistance during 10 Be measurements at the ASTER accelerator mass spectrometry (AMS) national facility (CEREGE, Aix-en-Provence, France). Radiocarbon dating was performed by the LMC14 Artemis AMS. This work was supported by the CNRS-INSU and CNES/ISIS research programs in France and by the Department of Geology and Mines in Bhutan.


Klaudia Bartolić Sirotić

Gaetano Benčić

Giulia Boetto

Anton Divić

Ida Koncani Uhač

Davor Munda

Marko Uhač

Sewn boat from Poreč, Istria (Croatia)

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

At the intersection of Obala Maršala Tita and Cardo Maximus streets, at a depth of -1,55 meters below the current mean sea level, a sewn boat was discovered lying on the ancient sea bottom together with some interesting harbor structures.

During April, excavation, in situ documentation and study, and salvage of the boat were carried out in collaboration between the Museum of the Poreč territory, the Archaeological Museum of Istria, the Directorate for Protection of Cultural Heritage of Croatia and the Centre Camille Jullian (Aix Marseille University, CNRS). Companies Macuka Ltd. and Kapitel Ltd. provided technical support for the archaeological excavation. After desalinization, the boat was send few days ago to Laboratory Arc-Nucléart of Grenoble for the conservation process. Before, it will be the object of dendrochronological analyses and detailed documentation.

The boat (named Poreč 1) was found adjacent to the remains of the Roman-era waterfront made of large stone blocks, following an approximate East-West orientation parallel to the current shoreline. The boat was covered with large stones. The western extremity was better preserved than the eastern part. Abundant ceramic and amphorae material were also found. Its preliminary study gives a date for the abandonment of the boat in the middle of the 1 st century AD. Its construction could be hypothetically placed in the first half of the 1 st century AD.

After fine cleaning, the boat was documented and prepared for various stages of 3D documentation by digital photogrammetry. Structure was also sampled to identify the species of wood used in construction.

The Poreč boat was preserved in the length of 5,.7 m and width of 1.4 m. Its longitudinal section preserved some of the original curvature in its western part and in the transversal section amidships.

The boat remains consisted of a keel and 2 transitional timbers, a total of 13 strakes (7 along the northern side and 6 along the southern side), 9 floor-timbers, one keelson/mast step, and 2 stringers.

Based on the position of the mast step and form of the mast recess it was possible to establish that the prow of the boat was on the eastern side. Unfortunately, this part of the boat was particularly damaged while the stern was in a better state of preservation.

The keel, 5.6 m long, was 5 cm sided and 6 cm moulded, and it had a rectangular, non-chamfered section. The hook scarfs on boat ends connected the keel with transitional timbers. The scarfs were secured by vertical treenails and a metal nail.

The planking was flush-laid and assembled by stitching. Three strakes were identified as stealers (PN2, PN3 and PS3), while three others presented repairs (PN4, PN5, PS2). Planks were 8 mm thick on average, with a maximum width for the northern garboard (PN1) of 15.5 cm (garboard PN1). The average with of the strakes was 11.7.

Planks were connected with simple over-edge stiches (parallel sewing pattern). The channels were perpendicularly drilled through the thickness of the planks. The outer internal edges of the planks bear traces of small rectangular recesses to protect the stiches.

The seams were blocked by small tapered pegs, with an average diameter of 7 mm. The pegs were driven from the inner side of the hull at an average distance of 2.6 cm. Some repair pegs covered by the floor timbers were driven from outside of the hull. A wadding pad was placed on top of the seams, and held tightly by the ligature. A thick layer of pitch completed the waterproofing system.

Nine floor timbers were preserved in place, and traces of additional ones were identified due to the rows of treenails still visible on the planking. Rectangular in section, on average 5 cm sided and 7 cm molded, they were spaced 38 cm apart, except between floor timbers F15 and F16, where the spacing was twice (80 cm). The floor timbers were connected to the planking by tapered treenails.

The lower face of the frames was crenellated in order to avoid crushing the wadding pad. On the bottom of the hull, the openings were quadrangular in shape and oversized (measuring approximately 6 cm in width and 2 cm in height), having the additional function of limber holes.

A keelson/mast-step timber, 1.3 m long, 11 cm wide and 5 cm thick, was installed on the frames between F10 and F13. It was found broken in two parts. This timber is wider around the mast recess, narrowing towards the extremities. Recesses were carved into the lower surface in order to be settled on frames.

The mast recess is rectangular in shape, 8 cm long, 4 cm wide and 4 cm deep. Its eastern side has been cut obliquely, at an angle of approximately 85º, in order to facilitate raising and lowering of the mast.

Finally, the two stringers were attached to the floor timbers by iron nails.

The Poreč boat shares some characteristics with others sewn boats found in Istria (Pula 2) and Dalmatia in Caska (Caska 1, 3 and 4) and Zaton (Zaton 1 and 3). The structure, sewing pattern and waterproofing system is similar. Preliminary examination of the wood of the Poreč boat testify the use of deciduous species of wood, in particular beech of planking as attested in Caska and Pula.

The Poreč boat was built according to a shell principle, and its shape was based on a longitudinal strake-oriented principle. Methods were shell-first.

The boat, long in origin around 7 m, was propelled using both sail and oars, and was used for coastal transportation of goods, peoples and animals.

In conclusion, since the first discoveries of sewn ships in Adriatic, new and significant ship-finds have improved the corpus offering new and valuable data to define the local traditions in shipbuilding that, originated in pre-Roman Istria and Liburnia, survived after the Roman conquest in this part of the Northern-eastern Adriatic nautical space.

Recent discoveries from Istria and Dalmatia strengthen the hypothesis of the survival during the Roman era of various traditions of sewn ships in the Adriatic region -in North-western Adriatic as well as in the North-eastern Adriatic. These traditions are marked not only by the construction of vessels structurally and functionally different, but also by the use of different sewing patterns or mixed techniques, the choice of different types of wood and raw materials for the hull, the assembly and the waterproofing system.

The Adriatic tradition of sewn ships has a local origin and it is strongly connected to the native pre-Roman shipbuilding tradition, as confirmed by the discovery of the Zambratija sewn boat dating back to the end of the Bronze age.

1

  Text of the communication "Sewn boat from Poreč, Istria (Croatia)" presented at the Sixteenth International Symposium on Boat & Ship Archaeology (ISBSA 16) Zadar, Croatia 26 September -1 October 2021 by Klaudia Bartolić Sirotić, Museum of the Poreč Territory, Poreč (Croatia) Gaetano Benčić, Museum of the Poreč Territory, Poreč (Croatia) Giulia Boetto, Aix-Marseille University, CNRS, Centre Camille Jullian, Aix-en-Provence (France) Anton Divić, Aix-Marseille University, CNRS, Centre Camille Jullian, Aix-en-Provence (France), Ida Koncani Uhač, Archaeological Museum of Istria, Pula (Croatia) Davor Munda, Museum of the Poreč Territory, Poreč (Croatia) Marko Uhač, Ministry of Culture, Directorate for Protection of Cultural Heritage, Conservation Department in Pula, Pula (Croatia) In the Adriatic region numerous examples of sewn boats have already been documented. This map shows an already rich archaeological corpus of sewn boats recently supplemented thanks to archaeological excavations in the antique colony of Parentium. At the beginning of 2020, the municipality of Poreč started to renovate the southern city waterfront. The archaeologists of the Museum of the Poreč territory were in charge of the archaeological supervision and excavation.


Anton Divić

Giulia Boetto

Krunoslav Zubčić

Zubčić Krunoslav

Digital reconstruction hypothesis of the Roman-era fluvial barge from Kamensko, Croatia

Text of the communication "Digital reconstruction hypothesis of the Roman-era fluvial barge from Kamensko, Croatia" presented at the Sixteenth International Symposium on Boat & Ship Archaeology (ISBSA 16) Zadar, Croatia 26 September -1 October 2021 by Anton Divić, Aix-Marseille University, CNRS, Centre Camille Jullian, Aix-en-Provence (France), Giulia Boetto, Aix-Marseille University, CNRS, Centre Camille Jullian, Aix-en-Provence (France) Krunoslav Zubčić, Hrvatski restauratorski zavod / Croatian Conservation Institute, Zagreb (Croatia) This presentation is focused on the reconstruction process of the Roman flat-bottomed ship found in the Kupa River in Kamensko. After a synthetic overview of the main architectural characteristics of the barge as already presented in the fifteenth edition of the ISBSA held in Marseille in 2018, we will address the theoretical principles behind the hypothesis of reconstruction of the shape and structure, explaining why and how some choices were made. Following will be the presentation of the digital techniques and methods used and, finally, an overview of the hydrostatic data extracted from the 3D model.

Kamensko shipwreck is that of a flat-bottomed Roman-era fluvial barge dated between the 1st century BC and 3rd century AD. Its remains, 2 m wide and preserved in the length of 12 m, are lying on the riverbed of the Kupa in Croatia, together with its relatively intact cargo of several hundred bricks of the type "Lydion". The shipwreck presents interesting structural features that place it in close relation with several other Roman-era shipwrecks from the Danube hydrographical basin.

The preserved remains of the barge consist of two massive hewn bilge strakes, between which 2 bottom strakes made of 4 planks have been inserted. A forwards-placed rectangular longitudinal mast step has been preserved in situ, as well as 5 pairs of frames attached to the strakes by treenails. Additional joining was achieved with tightly-spaced metal clamps, driven both from the outside and the inside of the hull.

The shipwreck was discovered in 2009 and was systematically excavated between 2015 and 2018 by a team from the Croatian Conservation Institute and the Centre Camille Jullian. Each step of the excavation was followed by meticulous documentation in the form of numerous measurements, drawings, detailed photographs, generation of photo-mosaics and three-dimensional photogrammetric models, continued during the post-excavation phase through the creation of both two-dimensional and three-dimensional site plans.

It was this abundance of high-quality documentation, in combination with the high level of preservation of the barge's structure, with approximately 54% of the original structure preserved, that allowed its hull to be reconstructed with a fair level of certainty.

The first phase in the reconstruction process was the formation of a theoretical reconstruction hypothesis. This theory of what the barge must have looked like during its last voyage and how it was navigated was not conjured in one instance, but was being developed constantly throughout the four excavation campaigns and numerous hour of post-excavation documentation and analyses, even if remaining only in a theoretical form. Besides through the study of the documented remains of the shipwreck, the original form of the barge was being reconstructed through comparison with archaeological, iconographic, as well as ethnological parallels.

One of the greatest challenges was the attempt to reconstruct the original length of the barge. Two methods were applied, both of them resorting to analogies with preserved Gallo-Roman fluvial vessels.

The first method for reconstructing the original length of the barge took into account the beam-tolength ratio of several well preserved Roman-era fluvial barges. This method was made possible due to the fact that the beam of the Kamensko barge that was conserved in its entirety. For purposes of comparison, a total of 11 Gallo-Roman fluvial vessels were taken in consideration, whose state of conservation and availability of publications allowed for a viable reconstruction of their original size. This resulted in an average beam-to-length ratio of 1:7.4. When multiplied with the 2.10 m beam of the Kamensko barge, its original length can be reconstructed to 15.5 m.

Second calculation introduced in order to validate, deny, or modify the obtained length took into account the position of the barge's mast step in relation to its prow, both of which were preserved in good condition on their original positions. The towing mast on fluvial barges is generally placed well forwards, approximately at 1/3 length from the prow, in order to keep the vessel parallel to the riverbanks and facilitate steering while being towed. The distance between the mast step mortise and the prow in relation with the full length of the vessel was calculated on the example of three Gallo-Roman fluvial barges, all resulting in an almost identical ratio. This ratio, applied in the case of the Kamensko barge, allowed us to propose a second hypothesis of its original length of 16.5 m.

If we decide to apply the middle value of the two close calculations, the original length of the Kamensko barge can be reconstructed around 16 meters.

The missing stern of the vessel was reconstructed as a mirrored version of the well-preserved prow. This decision was influenced by the symmetry that is often found in the number of Gallo-Roman vessels, as well as by the lack of possibility to introduce any kind of evidence that would point otherwise.

As it is very often the case with ancient shipwrecks, both the propulsion and direction system of the Kamensko barge did not survive the wrecking and site formation processes and will need to be reconstructed through deduction as well as through comparison and analogies with their counterparts documented in other archaeological, iconographical and ethnographical sources.

The only surviving element of the Kamensko barge connected to its propulsion is the mast step. However, the first question is what kind of mast did it support: a towing mast, a sailing rig, or both? In short, the possibility of the sail on the Kupa was dismissed due the several reasons: the relatively small width of the river, numerous developed meanders; absence of prevailing winds; no use of sail documented in the ethnographic record. On the contrary, this record testifies the large use of hauling.

These constrains, together with some archaeological, iconographic and textual sources that bear witness to the use of towing as a propulsion method employed on European rivers during the Roman period, allowed for a conclusion that the mast housed in this mast step was a towing one.

The evidence that allowed us for a more detailed reconstruction of this propulsion system came in form of three rectangular notches impressed into the upper face of the mast step, interpreted as traces of three planks that constituted the ensemble supporting the mast in place. Those would require a transversal support in form of a beam, possibly similar to the one recoded at the contemporary Zwammerdam 4 shipwreck, whose secondary function would be to provide the barge with additional transversal support.

For the reconstruction of the missing directional system, the archaeological, iconographical, and ethnographical parallels were once again introduced, speaking in favor of an axially positioned steering oar, complemented with a stern platform that would facilitate its handling and provide the steersman with a better line of sight. The eventual use of oars to keep the vessel in the channel when sailing downstream is also possible.

With the theoretical reconstruction hypothesis established, next stage in the Kamensko barge reconstruction process consisted of its materialization. In accordance with the latest approaches that are becoming standardized within the framework of nautical archaeology, an exclusively digital approach was undertaken, resulting in a digital, although functional and full-size restitution that can be tried and tested in a repeatable fashion. The barge's architectural structure and technical system were therefore three-dimensionally modelled in Rhinoceros 3D software.

First step in this process was the creation of a 3D model of the preserved remains. This is one of important steps of the process as the archaeological remains will constrain all the following steps.. This was achieved mostly with the help of georeferenced three-dimensional photogrammetric models of the shipwreck that served as a base layer for the 3D drawing of the remains in the Rhinoceros 3D software. The photogrammetric models are made out of a dense triangle polygon meshes, with the triangle points serving as attachment points for the three-dimensional drawing. In this way, the linear drawing of the shipwreck was effectuated that was to serve as a starting point for the subsequent modelling of surfaces, that is, for the rendering of the shipwreck's structure.

Information that could not have been extracted directly from the photogrammetric models was supplemented with numerous measurements and sketches hand-taken during the excavations that proved invaluable to the creation of the three-dimensional remains' model and the subsequent reconstruction.

The creation of the three-dimensional remains' model includes also the parts of the cargo and ship structure that were not visible on the 3D model obtained from the photogrammetry, as the outer surface of the barge. The modeling of the hull's outer shell was effectuated according to the direct recording underwater (ex. Thickness) and on the pieces extracted from the riverbed (sampling of the hull). It also allowed us the extraction of preliminary information concerning the barge's tonnage and displacement. But most importantly, it functioned as a base on which the subsequent reconstruction was modelled, serving as a foundation on which the laid-out theoretical ideas were modelled into a physical form.

The starting point for the reconstruction was the amidships cross-section of the shipwreck that was deemed not to present any post-deposition deformities. The upper edge of the side of the barge's caprail of the well-preserved starboard side was raised for several centimeters due to the traces of wear. The flank of the port side, preserved up the turn of the bilge, was therefore modeled on its starboard counterpart in order to obtain the reconstruction of the amidships cross-section. This amidships cross-section, once established and modelled, served as a reference point according to which the abaft starboard flank, presenting significant levels of distortion, was remodeled. In practical terms, this was achieved by establishing multiple cross sections of the modelled elements and straightening them one by one back to their original position in accordance to the established crosssection. This method was applied for the entirety of the deformed starboard flank. With the starboard flank in place, the shape of its port counterpart was modeled according to it, following however the different dispositions of bottom planks. This exclusively digital reconstruction method, besides being relatively time and cost effective, allowed for the reconstruction of a full-sized vessel that can subsequently be used for simulating its cargo capacity, displacements, and navigational capabilities within its area of navigation. Every architectural element was modelled in real size and had its respective material and characteristics assigned to it, which allows us a straightforward calculation of their volume and weight. This, in turn, leads to the relatively simple calculation of the vessel's displacements and tonnage, while the combination of the obtained data and the reconstructed lines of the hull allowed precise hydrostatic analyses to be generated in naval architecture software, in this particular case Orca 3D. The hydrostatic studies, however simplified in comparison with those of seagoing vessels, allowed some interesting information to be extracted.

With its reconstructed dimensions of 16.3 m in length, 2.1 m in width, and 0.7 m in depth, the 3D reconstruction allowed the calculation of the volume of the barge's structure including the assembly elements to 2.45 m3. This, in turn, allowed the barge's light displacement to be calculated to 2 tons, presenting a navigation draft of 17 cm when navigating empty.

Although the majority of the barge's brick cargo remained preserved in situ, the abaft cargo disposition was unknown due to missing stern of the barge, while the reconstructed form and volume of the hull allowed the full cargo disposition to be reconstructed as well. As correctly presumed in the early stages of excavation, the barge was transporting a cargo of 600 bricks occupying a volume of 4.65 m3. With each brick weighting 11 kg on average, plus the 140 kilos of dunnage in form of tree trunks and branches, the barge's deadweight tonnage at the moment of the wrecking was calculated at 6.75 tons.

The combination of barge's deadweight tonnage and light displacement allowed for its loaded displacement to be calculated. During its last voyage, the Kamensko barge is presumed to have had a displacement of 8.84 tons, that would result in a navigational draft of 47 cm, exactly 2/3 of the hull's total height. Calculations aimed at determining the maximal displacement of the vessel were undertaken as well, by simulating conditions that would lead to gradual sinking of the barge of 5 cm in turn, while still maintaining a minimal necessary freeboard for safe navigation. Which, according to a comparison with a modern document showing sand dredgers from the Haute-Garonne, can be kept at a very minimum. However, if we play it safe and allow for 15 cm of freeboard, the maximal displacement of the Kamensko barge amounts to just over 13 tons, allowing it to transport 11 tons of deadweight tonnage.


Reza Razi
email: reza.razi@grenoble-inp.fr

Majid Mehrasa
email: majid.mehrasa@grenoble-inp.fr

Khaled Hajar
email: khaled.hajar@grenoble-inp.fr

Mehrdad Gholami
email: mehrdad.gholami@grenoble-inp.fr

Ahmad Hably
email: ahmad.hably@grenoble-inp.fr

Seddik Bacha
email: seddick.bacha@grenoble-inp.fr

Antoine Labonne
email: antoine.labonne@grenoble-inp.fr

Predictive smart charging of plug-in electric vehicles for parking fleet scenario with modular converters

In recent years, smart charging of electric vehicles has been considered in various scenarios for all parties involved. However, the issue of smart charging is less explored in the parking lot scenario than in other ones, where the number of vehicles is typically greater than the number of charging points. In addition, the charging period is temporary, and modular converters are used to reduce costs in this scenario. In this paper, a modified algorithm is presented to connect each vehicle to the suitable and pursue the desired goals. The effectiveness of the proposed method is evaluated through simulation.

INTRODUCTION

Nowadays, the environmental problems associated with traditional vehicles with fossil fuels are clear enough to everyone. In this regard, electric vehicles (EVs), as a new generation and replacement of old vehicles, also have challenges [START_REF] Husain | Electric Drive Technology Trends, Challenges, and Opportunities for Future Electric Vehicles[END_REF]- [START_REF] Al-Hanahi | Charging Infrastructure for Commercial Electric Vehicles: Challenges and Future Works[END_REF]. One of the main concerns related to these vehicles is the significant rise in energy demand [START_REF]Electric-Vehicle Smart Charging Innovation Landscape Brief[END_REF]. In fact, a similar pattern of users' behavior, the design of ultra-fast chargers, and the significant growth of plug-in electric vehicles (PEVs) will impose heavy pressure on power distribution networks. In other words, PEVs increase the energy density of the load profile at certain times, especially around early morning and early evening. Therefore, this is a fact that the lack of energy demand management of PEVs will become a serious challenge in the near future. From the point of view of the economy and electrical sciences, this congestion leads to higher prices, reduced efficiency, and increased investment costs. Therefore, in an ideal plan, all parties involved in the issue are trying to resolve it through cooperation. In this regard, smart charging of PEVs has become a hot trend in recent years [START_REF] Ovalle | Grid Optimal Integration of Electric Vehicles: Examples with Matlab Implementation[END_REF], [START_REF] Hajar | Smart charging impact on electric vehicles in presence of photovoltaics[END_REF]. In other words, various methods control the output power of PEV as a flexible element to optimize the desired objective functions. In a simple view, PEV can be charged when the price or load is low or can be discharged in the opposite situation when the price is high. In general, there are two main control methods for smart charging of PEVs including the rule-based method [START_REF] Turker | Rule-based charging of plug-in electric vehicles (PEVs): Impacts on the aging rate of low-voltage transformers[END_REF]- [START_REF] Zhou | Model Prediction and Rule Based Energy Management Strategy for a Plug-in Hybrid Electric Vehicle With Hybrid Energy Storage System[END_REF] and optimization method [START_REF] Turker | Optimal Minimization of Plug-In Electric Vehicle Charging Cost With Vehicleto-Home and Vehicle-to-Grid Concepts[END_REF]- [START_REF] Franco | A Mixed-Integer Linear Programming Model for the Electric Vehicle Charging Coordination Problem in Unbalanced Electrical Distribution Systems[END_REF]. It can be claimed that in recent research, optimization algorithms have received more attention due to the better results, which have been improved in realtime with the help of the second method [START_REF] Guo | Two-Stage Economic Operation of Microgrid-Like Electric Vehicle Parking Deck[END_REF]. Also, researchers try to employ model predictive control (MPC) to deal with uncertainties in the desired structure. In one of our recent works in [START_REF] Razi | A userfriendly smart charging algorithm based on energyawareness for different PEV parking scenarios[END_REF], user-friendly smart charging was proposed by considering the charging cost and battery health, in such a way that the user plays a key role in choosing the desired charging scenario for PEV. In another work in [START_REF] Mehrasa | Power Management of a Smart Vehicleto-Grid (V2G) System Using Fuzzy Logic Approach[END_REF], a Mamdani-based fuzzy strategy was developed to manage the output power of PEVs. In addition, a simple technique was presented in [START_REF] Razi | Limiting discharge cycles numbers for plug-in electric vehicles in bidirectional smart charging algorithm[END_REF] to limit the number of discharge cycles for the bidirectional smart charging scenario. However, the main algorithm was only used for one vehicle in residential or workplace applications. Therefore, in most of past literature, less attention has been paid to energy management of parking lot scenarios where the number of charging points (CPs) is usually less than the number of PEVs and the parking time is temporary. In this paper, a user-friendly predictive algorithm for smart charging of PEVs in the parking lot application is presented, which is based on the following salient features:  Developing a novel algorithm to assign each PEV to the appropriate CP  Compatibility with modular converters  Independent of the number of PEVs with plug-andplay capability  Minimizing the charging cost by maintaining the priority of PEVs and satisfying the restrictions  Interacting with the user and providing possible results related to the cost and state of health (SOH)

Paper n° 371

CIRED 2022 Workshop 2/5

In the rest of the paper, a typical structure of a parking lot charging station is described in the next section. After that, the proposed algorithm will present. Finally, the simulation results will be shown to verify the effectiveness of the proposed algorithm.

PARKING CHARGING STATION STRUCTURE

In this section, it is better to start with a general structure of the electrical elements in the parking lot charging station. A Typical structure of parking lot charging station is shown in Fig. 1. As it can be seen in this structure, there are several CPs, a utility grid, local loads, and photovoltaic (PV) arrays. In smart charging algorithms, an energy management system (EMS) unit communicates with all the elements. Power directions for different components show positive directions of their power flow.

The main application of this structure is in shopping malls, city centers, government offices, hospitals, restaurants, airports and universities. In this regard, many PEVs will come to the desired parking lot, which the number of PEVs may be more than the capacity of CPs. Usually the parking time and energy demand of vehicles in this application is limited. The number of each PEV only indicates the priority of arrival of that vehicle. Also, each PEV can be connected to an arbitrary CP where it is not possible to change that after specifying. Chandra Mouli et al. [START_REF] Mouli | Comparison of system architecture and converter topology for a solar powered electric vehicle charging station[END_REF] tried to connect several vehicles to a single CP simultaneously by presenting a modular converter topology. The desired topology is shown in Fig. 2. In fact, due to the high energy density of CPs and reduced investment costs, this topology was very successful. It is assumed thast in this structure, there are two CPs, which can charge three PEVs simultaneously while maintaining the maximum power of CP (Pmax-CP = 10 kW). However, energy optimization between different elements and CPs in the presence of this topology requires a predictive algorithm to assign each vehicle to the appropriate CP. At the same time, the number of PEVs may still be more than the input of modular converters. In this case, the proposed algorithm should be able to reserve and supply energy related to vehicles in the charging queue, which will be presented in the next section.

PROPOSED SMART CHARGING ALGORITHM

As mentioned earlier, the main idea of this paper is inspired by the previous work of the authors in [START_REF] Razi | A userfriendly smart charging algorithm based on energyawareness for different PEV parking scenarios[END_REF], which could not be generalized for parking lot scenarios. In this regard, the proposed smart charging algorithm is improved and its flowchart is presented in Fig. 3.

Initially, data such as electricity price, weather conditions, and the number of CPs is obtained in the first step. Then, the desired information is estimated in the second step or is received directly from the user in the third step. In fact, the PV generation power and load profile should be estimated and arrival time, departure time, the initial state of charge (SOC), and final SOC can be received. Although obtaining some information from the driver may not be very accurate, it greatly reduces uncertainty. After that, as it can be seen, the algorithm uses two scheduling methods, including predictive linear optimization in the fourth step and a real-time scheme in the seventh step. In parallel with the optimization in the fourth step, the appropriate CP is assigned to each vehicle in the fifth step. Then in step six, the user or driver will select one of the possible scenarios. Finally, the real-time controller is used to compensate for errors caused by uncertainties. However, in this paper, the seventh step is ignored because there is no real-time implementation. The following is a more complete description of linear optimization and the algorithm for assigning PEVs to CPs.

Here, unlike the previous work, all PEVs are controlled based on a single optimization controller. The proposed optimization is employed in the form of MPC, which its details are shown in Fig.

                    (1)
where n, i, Cgrid, Δt, EEV, λenergy and λbattery represent the total number of PEVs, the index of each PEV, the electricity price, optimization time step, the PEV battery energy, charging priority coefficient and battery efficiency coefficient, respectively. As it is clear, the second term of the objective function is related to the priority of charging.

In fact, a PEV has a higher energy coefficient when it arrives to the charging station earlier. This factor is defined in this work as follow:

, ( 1) 100 09; 1,...,

energy i n i E i n        (2)
Besides, the third term avoids useless charging and discharging of batteries due to the loss of energy as follow: , , , 1 ; 0 ; 1,..., 0.9 ; 0

EV i battery i EV i P i n P           (3)
Also, the power balance equation is:

On the other hand, the proposed scheme for assigning PEVs to CPs is working as shown in Fig. 5.

As it can be seen, this algorithm is repeated for all n vehicles to find their proper CP. The algorithm is activated only once for each vehicle when the EV enters the charging station (t = tarrival(i)). If the vehicle number i be equal to or less than the number of the CPs (NCP), the vehicle is connected to the corresponding CP (CPi = i). In this way, one-to-one correspondence is created between the EVs and the CPs, and the maximum dispersion is respected.

Otherwise, when the EV number is more than the number of CPs, there is a different procedure for selecting the CP. In this regard, for all CPs and for the time period between the arrival and departure times of the desired vehicle (tdeparture(i)tarrival(i)), the available power of CP (PoCPj) and its cost (CoCPj) are calculated. This variable represents the remaining power capacity of the CP, which is equal to the maximum power minus the power of previous vehicles connected to that CP. max

j available j CP j EV i ij PoCP P P P         (5)
In the next step, the cost-to-power ratio (CoCPj/ PoCPj) of each CP is obtained such that the CP with the lowest value is selected as the appropriate CP of the vehicle. In fact, the algorithm considers the CP with more capacity and lower price for EVs. When each EV is added to a appropriate CP, the following constraint is activated:

max EV new EV i CP j ij P P P        (6)
At the end, this procedure is repeated for the next EV. As a result, the power of vehicles, grid power, and associated CPs will be obtained.

SIMULATION RESULTS

In this section, the simulation results of various cases in the presence of proposed method will be analyzed. It is assumed that the grid power limitation is 30 kW. In this regard, the results for a parking lot with two CPs are shown in Figs. 6789based on the provided information in Table 1.

Here, due to the short parking time and constant electricity price, the unidirectional scenario is only considered such that the PEV will act as a flexible load.

In the first case study, the first EV arrives at the parking lot at 8 am. The simulation results are shown in Fig. 6, which includes the power of different elements, charge status and cumulative cost. As it can be seen, EV1 is immediately connected to the first CP and it is charged with solar energy to reach the final SOC. It should be noted that during this period, the energy exchanged with the grid is null.

As shown in Fig. 7, in the second case study, the EV2 is connected to the CP2. Although the CPs of the first and second EVs are different, EV2 is charged after the first vehicle to use the renewable energy produced inside the parking lot. Here, due to being at the PV power peak, the EV2 is charged with maximum power and at a constant rate.

Then, in the third case, the simulation results with the addition of the third vehicle are shown in Fig. 8. It should be noted that EV3 arrives at the parking lot earlier than EV2, but the second car has reserved CP2 by announcing earlier. Therefore, according to the proposed algorithm, EV3 is connected to CP1. In this case, as in previous simulations, no energy is received from the grid to charge EVs.

Finally, the simulation results in the presence of four vehicles are shown in Fig. 9. In this case, there is no enough photovoltaic energy to charge all the cars. Therefore, the grid utility needs to be involved in the charging process of EVs. Here, due to the provided data and limitations, the fourth EV is also connected to the CP1, such that the changes in the output power of other EVs are also evident. Therefore, as expected, with the addition of other vehicles, the system constraints, especially the capacity of the CPs, are activated and the optimal conditions change.

Conclusion

Smart charging of EVs in the parking lots in the presence of modular converters can have technical and economic benefits for the driver, station owner and network operator. For this purpose, in this paper, a modified predictive intelligent charging algorithm was presented that determines the ideal charging point and charging process of different vehicles. In this regard, in the proposed objective function, several goals were pursued, including minimizing the cost of energy exchanged with the grid, prioritizing for charging of EVs and considering battery efficiency. In addition, to assign the suitable charging point, the available power at all charging points and electricity cost were used. Finally, it was verified with the help of simulation results that the proposed algorithm has the flexibility to operate in the desired application.

Fig. 3 .Fig. 2 .

 32 Fig. 3. Flowchart of the proposed smart charging algorithm.

Fig. 1 .

 1 Fig. 1. Typical structure of parking lot charging station.

Fig. 5 .Fig. 4 .

 54 Fig. 5. Proposed flowchart of the CPs assigning (Fifth step).Stop

Fig. 9 .Fig. 8 .

 98 Fig. 9. Simulation results for four PEVs and two CPs.

Fig. 7 .

 7 Fig. 7. Simulation results for two PEVs and two CPs.

Fig. 6 .

 6 Fig. 6. Simulatiosn results for a single PEV and two CPs.

Acknowledgments

This work has been supported by the aVEnir project of the PIA operated by ADEME.