Avinash Kumar Singh

Mohamed Adjel
email: adjelmohamed@gmail.com

Vincent Bonnet
email: vincent.bonnet@laas.fr

Robin Passama

Andrea Cherubini

A Framework for Recognizing Industrial Actions via Joint Angles

This paper proposes a novel framework for recognizing industrial actions, in the perspective of human-robot collaboration. Given a one second long measure of the human's motion, the framework can determine his/her action. The originality lies in the use of joint angles, instead of Cartesian coordinates. This design choice makes the framework sensor agnostic and invariant to affine transformations and to anthropometric differences. On AnDy dataset, we outperform the state of art classifier. Furthermore, we show that our framework is effective with limited training data, that it is subject independent, and that it is compatible with robotic realtime constraints. In terms of methodology, the framework is an original synergy of two antithetical schools of thought: modelbased and data-based algorithms. Indeed, it is the cascade of an inverse kinematics estimator compliant with the International Society of Biomechanics recommendations, followed by a deep learning architecture based on Bidirectional Long Short Term Memory. We believe our work may pave the way to successful and fast action recognition with standard depth cameras, embedded on moving collaborative robots.

I. INTRODUCTION

In recent years, industry is adopting more and more collaborative robots (cobots). This opens the door to new scenarios, where human and robot share the same workspace [START_REF] Hvilshøj | Autonomous industrial mobile manipulation (aimm): Past, present and future[END_REF]. For a robot to be effective in this context, there must be a tight choreography between the two (human and robot). Typically, for successful human-robot collaboration, it is crucial that the robot understands human actions. This will be helpful in numerous scenarios and tasks including perilous environments, heavy object lifting, tool pick and place, etc. Human-robot collaboration should not only increase productivity; it should also reduce musculoskeletal disorders [START_REF] Schneider | OSH in figures: Occupational safety and health in the transport sector -an overview[END_REF] induced by repetitive tasks performed in awkward postures. Therefore, online recognition of the human actions will also provide feedback on ergonomic risks [START_REF] Malaisé | Activity recognition for ergonomics assessment of industrial tasks with automatic feature selection[END_REF].

Human action recognition has been an active research domain since more than two decades [START_REF] Rodríguez-Moreno | Video activity recognition: State-of-the-art[END_REF], [START_REF] Zhang | Rgb-d-based action recognition datasets: A survey[END_REF]. Some researchers focus on recognizing human actions from RGB and RGB-D videos [START_REF] Karpathy | Large-scale video classification with convolutional neural networks[END_REF] [START_REF] Tran | Learning spatiotemporal features with 3d convolutional networks[END_REF], while others utilize the position of the human joints (skeleton) [START_REF] Du | Skeleton based action recognition with convolutional neural network[END_REF] [START_REF] Liu | NTU RGB+D 120: A large-scale benchmark for 3d human activity understanding[END_REF], obtained from motion capture systems [START_REF] Barnachon | Ongoing human action recognition with motion capture[END_REF] [START_REF] Ionescu | Human3.6m: Large scale datasets and predictive methods for 3d human sensing in natural environments[END_REF]. There are many labeled datasets online [START_REF] Rodríguez-Moreno | Video activity recognition: State-of-the-art[END_REF][5] [START_REF] Liu | NTU RGB+D 120: A large-scale benchmark for 3d human activity understanding[END_REF], to support or verify methodologies, particularly those based on machine learning. Yet, most of these datasets are limited to daily life actions, such as making coffee, reading the newspaper, using the telephone etc. Few focus on tasks which take place in industrial scenarios, such as pick/place a tool, screw, etc. As per our knowledge, there are only two industry-oriented motion capture datasets: InHard [START_REF] Dallel | Inhard -industrial human action recognition dataset in the context of industrial collaborative robotics[END_REF] and AnDY [START_REF] Maurice | Human movement and ergonomics: An industry-oriented dataset for collaborative robotics[END_REF]. In InHard, humans perform all actions while standing, although this is not always the case in industrial setups. AnDY instead, covers a broader variety of postures (e.g., walking, sitting). Therefore, we decided to use AnDY dataset to test our approach.

In terms of methods, most recent works on human action recognition rely on deep learning. Researchers have used a variety of neural networks such as Convolutional Neural Networks [START_REF] Wang | Actions ˜transformations[END_REF][15] [START_REF] Liu | Skepxels: Spatiotemporal image representation of human skeleton joints for action recognition[END_REF], Recurrent Neural Networks [START_REF] Zhang | On geometric features for skeleton-based action recognition using multilayer lstm networks[END_REF][18] [START_REF] Du | Hierarchical recurrent neural network for skeleton based action recognition[END_REF], Graph Neural Networks [START_REF] Shi | Skeletonbased action recognition with directed graph neural networks[END_REF][21] and combinations of the above [START_REF] Ullah | Action recognition in video sequences using deep bi-directional lstm with cnn features[END_REF]. All the cited studies have addressed the problem in the Cartesian space (XY Z coordinates of the human joint centers). In our work, we also opt for deep learning, but propose the use of joint angle space rather than Cartesian space.

In the context of Human-Robot Collaboration, two important characteristics are the low quality of data (since the sensor is often on-board the robot) and the need for fast processing. Specifically:

• sensor data can be noisy, incomplete or blurred (due to occlusions, limited field of view or robot motion), • data must be processed and classified (e.g., actions must be recognized) at a framerate that is high enough to control a robot, at least at 5-10 Hz. These issues make the action recognition problem much more challenging in Robotics than it is in the field of Computer Vision, where the papers [START_REF] Zhang | Rgb-d-based action recognition datasets: A survey[END_REF][23] [START_REF] Yasin | Keys for action: An efficient keyframe-based approach for 3d action recognition using a deep neural network[END_REF] present solutions, which do not account for low quality data nor for framerate requisites. Besides, they do not easily generalize to:

• viewpoint (orientation and translation),

• subject-specific anthropometry,

• sensors other than cameras. To address these aspects, one should train the classifier on datasets covering many cases, at a huge cost in terms of acquisition, labeling and training.

Our solution consists in mapping the 3D Cartesian coordinates of the human skeleton joint centers, which can be derived from diverse sources (RGB-D camera, stereophotogrammetric device, inertial suits, ...) to his/her joint angles. The joint angles representation has many advantages:

• it is sensor agnostic, i.e., one can apply it to any motion capture device measuring the human joint Cartesian coordinates; such device may be a dedicated motion capture sensor (e.g., optical, inertial, etc) or the combination of a sensor and processing algorithm (e.g., RGB-D image stream processed by openpose 1 ), • it is invariant to affine transformations (scaling, translations, rotations), • it generalizes well across the anthropometry of different subjects, • it is less sensitive to outliers, which can be processed via standard joint angle limits (from biomechanical tables). These advantages are crucial in the perspective of machine learning, since a joint-based classifier will require a drastically smaller amount of data than the state-of-art Cartesianbased classifiers, to achieve comparable generality. In a nutshell, the contributions of this paper are:

• on AnDY dataset, we show that the joint angle representation yields higher accuracy than the Cartesian representation and than the state of art classifier [START_REF] Malaisé | Activity recognition for ergonomics assessment of industrial tasks with automatic feature selection[END_REF]; • because of the advantages cited just above, our joint angles classifier is accurate across all AnDY actions, effective with limited data and person independent; • because of the lower dimension of the feature vector, training is much faster with joint angles than it is with Cartesian coordinates; • the classifier computation time is compatible with realtime constraints (even considering the addition of a vision-based skeleton tracker). Last but not least, our work is an original synergy of two schools of thought which tend to compete in the field of engineering, those of the model-based and data-based communities. Indeed, by complementing biomechanics (Sec. II) with deep learning (Sec. III), we outperform the state of art [START_REF] Malaisé | Activity recognition for ergonomics assessment of industrial tasks with automatic feature selection[END_REF] in the recognition of industrial actions (Sec. IV).

II. THE PIPELINE FOR JOINT ANGLE ESTIMATION

A motion capture system usually measures Cartesian poses of the body segments. To estimate joint angles from these measures, we developed a Capture System Software Library (CSSL), composed of an offline and an online phase. The offline phase consists in creating a biomechanical model of each subject based on a Unified Robot Description Format -URDF. At each new sample of time, the online phase updates the joint angles of the model by calculating inverse kinematics. Both phases are detailled below.

A. Biomechanical Model of the Human

The proposed 23 Degrees of freedom (DoF) biomechanical model is based on the recommendations of the International Society of Biomechanics [START_REF] Wu | Isb recommendations for standardization in the reporting of kinematic data[END_REF]- [START_REF] Wu | Isb recommendation on definitions of joint coordinate systems of various joints for the reporting of human joint motion -part ii: Shoulder, elbow, wrist and hand[END_REF]. This way, we can cope with the differences between scientific and clinical communities, and use -in future work -the model in the field of industrial ergonomics.

Figure 1 shows the model with the position of the Joint Center Positions (JCP) and the type of corresponding mechanical joints. The model base is attached to the pelvis segment. The pelvis is connected to each leg by a 3 DoF ball hip joint and the knee is represented by a hinge joint. The offline phase takes as input constant parameters the subject's gender, height, and mass and it estimates -using anthropometric tables [START_REF] Dumas | Estimation of the Body Segment Inertial Parameters for the Rigid Body Biomechanical Models Used in Motion Analysis[END_REF] -each segment's geometric and inertial parameters. With these, we generate the human model, as a Unified Robot Description Format -URDF. We do this with a C++ library (https://gite.lirmm.fr/ humar/humar-joints-estimator) with classes representing each part of the biomechanical model (segments, joints, limbs) and its parameters (mass, length, inertia and center of mass)2 .

B. Real-time Inverse Kinematics Estimation

Once the URDF has been generated, the pipeline's online phase uses at each iteration t the JCP, to estimate the corresponding joint angles of the subject. More formally, our inverse kinematics problem consists, at each iteration t, in estimating the joint angle vector θ t ∈ R 23 , corresponding to the measured 3D Cartesian coordinates of the 16 JCP, χt ∈ R 48 . The Forward Kinematics Model (F KM ), calculated from the human URDF, maps one representation to the other:

χ t = F KM (θ t ) . (1)
We can fit this model on the measures, by solving a constrained optimization problem, which consists in minimizing the distance between the measured and modeled JCP, while accounting for a regularization term:

θ * t = argmin θt∈R 23 ||F KM (θ t ) -χt || 2 2 + β||θ t -θ * t-1 || 2 2 subject to θ -≤ θ t ≤ θ + (2)
where β = 1 -3 is a weight allowing to avoid discontinuities, and θ + and θ -are the upper and lower joint limits, respectively. We solve (2) efficiently using the C++ library Ipopt [START_REF] Wächter | On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming[END_REF] along with CppAD library [START_REF] Bell | Algorithmic differentiation of implicit functions and optimal values[END_REF].

III. ACTION RECOGNITION

A. Dataset Description

To validate our method, we used AnDY dataset [START_REF] Maurice | Human movement and ergonomics: An industry-oriented dataset for collaborative robotics[END_REF], [START_REF] Maurice | Andydatalab-oneperson[END_REF]. Six industry-specific tasks (Screw High, Screw Middle, Screw Low, Untie Knot, Carry a 5kg load and Carry a 10kg load) are combined in 6 alternative sequences. Thirteen healthy adults (9 males and 4 females, of height 175.4 ± 7.9cm, weight 72.3 ± 14.4kg and age 25.7 ± 5.0years) participated in the study. Each participant repeated 5 times 3 sequences selected randomly among the 6. The dataset contains three types of human motion measurements: inertial motion capture data (Xsens MVN Link system, 240Hz, Xsens, Enschede, The Netherlands), stereophotogrammetric data obtained with a motion capture system (Qualisys, 120Hz, Goteborg, Sweden), and videos recorded by two RGB cameras [START_REF] Malaisé | Activity recognition for ergonomics assessment of industrial tasks with automatic feature selection[END_REF]. Within the 13 × 5 × 3 = 195 sequences of the dataset, the authors have labeled 8 actions: Reach, Pick, Place, Release, Carry, Fine Manipulation, Screw, Idle, defined in Table I.

Figure 2 shows the concatenation of actions over time, for one sequence of 107 seconds. Actions are color-coded and represented by circular sectors of radius proportional to their duration (indicated in seconds). For instance, the largest orange sector on the right indicates that the participant was idle at the end of the sequence (from 80.7 s to 107 s). We can see that some actions (Id, Fm, Re and Rl) are more frequent and on average longer than others (Ca, Pl, Pi and Sc). This is the case for all sequences, making the dataset distribution skewed, i.e., some classes have more samples than others.

B. From inertial suit measures to joint angles

This section describes how we utilized the AnDY dataset to build the training and testing subsets for our human action classifier. To this end, we exploited the data recorded with the inertial suit at 240Hz. At each iteration t, the inertial suit estimates the Cartesian Coordinates of 23 joints of the human body, denoted ψ t ∈ R 69 and listed in Table II. However, to calculate the corresponding joint angles with the CSSL, we require the 16 JCP indicated in Sec. II-B. To this end, we map ψ t to χ t , as shown in Table II. The head motions were not considered in our study, since they are prone to variations reflecting the subject's cognitive state. The so-called shoulder joint provided by the inertial suit corresponds to the center of the sterno-clavicular joints, and can be estimated only using regression methods. Therefore, we only use the cervical joint as center of rotation of the clavicules, as recommended in biomechanics [START_REF] Wu | Isb recommendation on definitions of joint coordinate systems of various joints for the reporting of human joint motion -part ii: Shoulder, elbow, wrist and hand[END_REF]. We do not consider the positions of the toes, and thus the ankle joint angles, and the Cartesian positions of vertebrae T12 and T3, since they displayed small amplitude in the investigated tasks and are therefore error prone. At each sample time, we set the model base at the pelvis, with the transformation matrix calculated from the upper legs, pelvis and L5 Joint Coordinates, measured by the inertial suit [START_REF] Cappozzo | Position and orientation in space of bones during movement: Anatomical frame definition and determination[END_REF]. Figure 3 shows the results of the CSSL for four representative postures of the investigated tasks.

C. Architecture of the Action Classifier

In this section, we describe the architecture of the classifier, which we have designed to map a sequence X of T consecutive feature vectors: to one of the 8 action classes in Table I. The size n of feature vector can be either 69 -when using all the JCP provided by the inertial suit, ψ -or 23 -when using the joint angles θ estimated using the CSSL. At each iteration k, our classifier assigns an action label to x k ; after T iterations, the classifier labels the sequence X t with the class which has obtained the majority of labels. Our approach is inspired by [START_REF] Zhang | On geometric features for skeleton-based action recognition using multilayer lstm networks[END_REF], [START_REF] Du | Hierarchical recurrent neural network for skeleton based action recognition[END_REF]. The authors used stacked LSTM for action recognition, but in this work, we tried [START_REF] Schuster | Bidirectional recurrent neural networks[END_REF], a bidirectional variant of LSTM (Bi-LSTM). Bi-LSTM adds one more LSTM layer, which reverses the We believe Bi-LSTM is particularly useful when dealing with action sequences, since the classifier can refer to future actions to label the current one. Figure 4 illustrates our Bi-LSTM architecture. We stacked three Bi-LSTM layers, to make the network preserve the spatial relation between joints. Each of the 3 hidden layers is composed of 32 cells (h = 32). Before training and validation, input x t is normalized (shifting it by its mean, and then dividing by its standard deviation). Then, the first layer applies to this normalized input. The output of the LSTM layer is then passed to the following layers. To avoid overfitting, we apply a dropout of 0.5 to the third layer. The stacked layers are followed by the final fully connected layer, with softmax activation to map x t to one of the 8 classes. We used categorical cross entropy as the loss function, Adam as optimizer, and we set the learning rate to 0.001.

X t = {x t-T , . . . , x t-1 } x k = {ψ k , θ k } ∈ R n (3)

IV. VALIDATING THE ACTION CLASSIFIER

In this Section, we present the tests that we ran to validate the classifier, developed with TensorFlow 2.7, on a computer with Intel © Core i7-1165G7@2.80GHz, 16 GB RAM and 4 GB GPU. The python code implementing the classifier can be found at https://gite.lirmm.fr/humar/ research-projects/humar_dnn. A video presenting the results is attached to this paper and available at https: //youtu.be/2GNWYsOfDYw.

As in [START_REF] Malaisé | Activity recognition for ergonomics assessment of industrial tasks with automatic feature selection[END_REF], we use a time window of 1 second with 50% overlapping, to define each sample of an action. Since data is measured at 240 Hz, this corresponds to T = 240 in [START_REF] Malaisé | Activity recognition for ergonomics assessment of industrial tasks with automatic feature selection[END_REF]. At this stage, some classes have many more samples than others, due to their frequency and execution length (Sec. III-A). Hence, to balance the data, we randomly discard some samples from the most represented classes (Id, Fm, Re, Rl).

In a first experiment, we randomly split the dataset into 70% training and 30% validation subsets, and validate it with the action classifier presented in Sec. III-C. We use the same architecture and hyperparameters for both representations of the input data, ψ ∈ R 69 and θ ∈ R 23 . Table III compares the F1 scores of our classifier with the two inputs ψ and θ, with the state of art Hidden Markov Model from [START_REF] Malaisé | Activity recognition for ergonomics assessment of industrial tasks with automatic feature selection[END_REF]. Both the figure and table show that the joint angle representation (θ) slightly prevails over the Cartesian one (ψ). Despite the small difference in accuracy, it should be noted that training is much faster with θ than it is with ψ. On average, a training epoch takes 130.13 s with θ and 230.38 s with ψ. This is because of their different dimensions (23 vs 69). In terms of accuracy, both θ and ψ outperform [START_REF] Malaisé | Activity recognition for ergonomics assessment of industrial tasks with automatic feature selection[END_REF], where the F1 score was 86.63% (see Table III). We believe that the reason is that Recurrent Neural Networks can handle long term data dependencies better than Hidden Markov Models, where the state at time t is related only to the state at t -1. Another advantage of our approach over [START_REF] Malaisé | Activity recognition for ergonomics assessment of industrial tasks with automatic feature selection[END_REF], is that it requires solely the joint positions, which can be measured with a variety of human motion capture devices. Instead, classifier [START_REF] Malaisé | Activity recognition for ergonomics assessment of industrial tasks with automatic feature selection[END_REF] needs 11 features, including higher-order kinematic variables (e.g., angular velocities and accelerations), which cannot be measured by most devices (e.g., depth cameras). Figure 5 is the confusion matrix of our Bi-LSTM with joint angles θ as input. For each action the values on the diagonal are the percentage of true positives, while the others are the percentage of misclassifications (false negatives or false positives). The very high values on the diagonal confirm the quality of our classifier across all eight actions.

To verify the performance across different subjects, we tried leave-one-out cross validation for joint angles input θ. We excluded all 15 trials of one participant from the training set, leaving data of the other participants. Then, we validated only with the excluded participant. We show the results for each action in Fig. 6, with the excluded participants in abscissa, and the F1 score in ordinate. For example, the first column indicates the F1 score of participant p5124 with the classifier trained on the rest of the population. Validation accuracy is similar (and high) across all actions for all participants, except p3327 (column 9). The reason could be that participant p3327 is left handed, while most participants (11 out of 13) are right handed. Ideally, the dataset should be equally distributed among left and right handed participants. Table IV shows the F1 scores when validating with all participants, the classifier was trained with and without participant p3327. For all actions, the F1 score is greater than 89.8%, and when excluding p3327 its decreases by 2.0% (average across the 8 actions). These results show the robustness of our classifier, with regards to new data. In a more challenging experiment, we excluded almost half (6/13) of the population and we trained the classifier with the rest. We excluded all samples of 4 male participants (p909, Fig. 6: Leave-one-out F1 score of our Bi-LSTM with input θ. Abscissa: excluded participants, ordinate: F1 score with that participant. p3327, p5124 and p9266) and 2 female participants (p8524 and p5535) from the training, while we used only these 6 participants for validation. We repeated the experiment with both inputs: raw inertial suit JCP ψ and joint angles θ. We report the F1 scores for all eight actions in Table V. Obviously, with less examples for training, all values have decreased with regards to Tables III. Yet, there is a clear gain in using joint angles θ compared to ψ. The average F1 scores were 85.48% and 68.19% for θ and for ψ, respectively. This shows that joint angles classification is robust to the use of limited data, and that it is subject independent. It also means that when using a biomechanical model and calculating joint angles, we make the data agnostic to the subject's biometrics. Let us conclude on the possibility of using our classifier for real time robot control in the context of human-robot interaction. This paper introduces the foundations of our future work, which aims at fast action recognition with low cost sensors embedded on a -possibly moving -robot. Ultimately, we believe these will be RGB-D cameras with the inclusion, in the perception pipeline, of a skeleton tracking algorithm such as Openpose. Such algorithms nowadays require long computing time (e.g. 60 to 80ms from our experience) even with the GPU embedded on our robot (NVIDIA GeForce GTX 1080). This will substantially reduce the data frame rate, from the current 240Hz, to approximately 10Hz. To verify the feasibility at such framerate, we subsampled the inertial suit data from the original 240Hz to the targeted 10Hz. This corresponds to taking a measure every 24 and setting T = 10 in (3). Under such conditions, we train and validate the classifier with both ψ and θ, to obtain the F1 scores shown in Table VI, along with the ones at 240Hz (copied from Table III). The results are very encouraging, with a loss of only 3% when subsampling the joint angles θ. Along with the accuracy, one should verify the feasibility in terms of computation time. Fulfilling the 10Hz constraint in the worse case scenario -camera acquisition and skeleton tracking taking 80ms -leaves 20ms for the other operations. These are: classification (i.e., the time required by the Bi-LSTM to label a sequence) and (when using joint angles θ) inverse kinematics estimation, described in Sec. II-B. Table VI shows the computation times of these two operations for the four versions of the classifier. These are averages measured on the mentioned computer, for one subject. The results show that even with the addition of inverse kinematics, the pipeline should not violate the 100ms constraint.

V. CONCLUSIONS

We have introduced and validated a framework for recognizing industrial actions from human motion data. We argue that relying on joint angles rather than on the Cartesian coordinates -commonly used in the literature -enhances the classification performance at many levels. The classifier is less sensitive to subject variance, faster to train, and compatible with robotic real-time constraints. Our claim is confirmed by the results obtained on inertial suit measures from AnDy dataset, to classify 8 different actions.

In the future, we will apply the classifier to data from low-cost depth cameras, combined with a visual skeleton tracker (e.g., openpose). Currently, to generate the biomechanic model of each subject, we must know his/her gender, height and weight. We will try to avoid this, by relying solely on sensor data. We also believe that the joint angles representation will help address challenges present in vision systems, such as occlusions and undetected joints.

1

  https://github.com/CMU-Perceptual-Computing-Lab/ openpose

Fig. 1 :

 1 Fig. 1: Representation of the 23 DoF human biomechanical model and relevant JCP. Estimated centers of masses are represented with red spheres.

Fig. 2 :

 2 Fig. 2: Concatenation of actions for one sequence in AnDY (Participant 541 Setup A Seq 3 Trial 1, see[START_REF] Maurice | Andydatalab-oneperson[END_REF]), starting from the mid left (Id 0.0-7.04). Each action is represented by a circular sector of radius proportional to its duration and labeled with the start and end time (in seconds). The color code is: orange (Idle), green (Reach), blue (Fine Manipulation), purple (Pick), black (Carry), magenta (Place), red (Screw), cyan (Release).

Fig. 3 :

 3 Fig. 3: Representative postures showing the results of the CSSL while: (a) standing in a neutral posture (Idle), (b) bending to pick-up an object with both hands (Pick), (c) walking while carrying an object with both hands (Carry), (d) bending to place an object with one hand (Place).

Fig. 4 :

 4 Fig. 4: Architecture of our action classifier. The input (blue) is n-dimensional vector x = {ψ, θ}, with n = 23 or 69. This is normalized, passed through 3 Bi-LSTM layers (orange), followed by a fully connected layer (green) with softmax activation to obtain the action class label (Y).

Fig. 5 :

 5 Fig. 5: Confusion Matrix of our Bi-LSTM with input θ.

TABLE I :

 I Definition of the actions from the AnDY dataset.

TABLE II

 II

: Mapping from the 23 joints (ψ) measured by the inertial suit to the 16 JCP (χ) needed by the CSSL for inverse kinematics.

TABLE III :

 III

F1 scores and input dimensions (n) of: our Bi-LSTM with input ψ, our Bi-LSTM with input θ and the state of art Hidden Markov Model

[START_REF] Malaisé | Activity recognition for ergonomics assessment of industrial tasks with automatic feature selection[END_REF]

.

TABLE IV :

 IV F1 scores (in %) of the θ classifier across all participants, trained with/without participant p3327.

TABLE V :

 V F1 scores (in %) of the classifier, trained with 7 of the 13 participants, with inputs ψ and θ.

	Performances	ψ@240Hz θ@240Hz ψ@10Hz θ@10Hz
	F1 score	98.99	99.42	95.44	96.38
	Classification time 1.70	1.44	0.18	0.17
	IKE time	N.A.	5.25	N.A.	5.25

TABLE VI :

 VI

F1 scores (in %) and average times (ms) needed for Classification and Inverse Kinematics Estimation (IKE), with four versions of the classifier.

Although in this work we focus on kinematics, the body segment inertial parameters encode information about the dynamics, which could be exploited in future work.

Pick up an object: starts when touching the object, ends when the arm stops moving with respect to the body.Place (Pl)Place an object: similar to Re, but with an object in hand.Release (Rl)Bring arm back after manipulation.Carry (Ca)Carry an object: starts at the end of Pi, ends at the beginning of Pl.Fine manipulation (Fm)Manipulate an object with dexterity.Screw (Sc)Rotate the hand with screwing motion (particular case of Fm). Idle (Id) Don't move the hands.

ACKNOWLEDGMENT

This research was carried out in the context of the SOPHIA project, which received funding from the EU Horizon 2020 research and innovation program under Grant Agreement No. 871237. We would like to thank the LARSEN team at INRIA for providing the AnDy dataset and Raphael Dumas from IFFSTAR-University of Lyon for his help on biomechanical modeling.


David Kaniewski
email: david.kaniewski@univ-tlse3.fr

Nick Marriner

Giovanni Sarti

Duccio Bertoni

Marco Marchesini

Veronica Rossi

Anna Lena

Alexandra Bivolaru

Majid Pourkerman

Matteo Vacchi

Rachid Cheddadi

Thierry Otto

Frédéric Luce

Daniela Cottica

Christophe Morhange

Northern Adriatic environmental changes since 500 AD reconstructed at Aquileia (Italy)

Keywords: Aquileia, Fluvial harbour, Ecosystems, Anthropogenic activities, Climate, Late Holocene, Italy

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

Introduction

The reconstruction of long-term environmental dynamics from ancient harbours has highlighted the importance of bio-and geo-sciences in understanding coastal man-made environments resulting from the Mediterranean's long and rich maritime history (e.g. [START_REF] Marriner | Geoscience of ancient Mediterranean harbours[END_REF][START_REF] Marriner | Ancient harbour infrastructure in the Levant: tracking the birth and rise of new forms of anthropogenic pressure[END_REF]. In Italy, numerous seaports and associated fluvial harbours have been intensively studied, providing an in-depth understanding of the human activities and anthropogenic pressures on environments since the emergence and development of these port complexes. Harbours such as Altinum-Venice [START_REF] Ninfo | The map of Altinum, ancestor of Venice[END_REF], Portus Lunae-Luni [START_REF] Bini | Palaeogeographies of the Magra Valley coastal plain to constrain the location of the Roman harbour of Luna (NW Italy)[END_REF], Portus-Rome [START_REF] Sadori | The plant landscape of the imperial harbour of Rome[END_REF][START_REF] Sadori | Archaeobotany in Italian ancient Roman harbours[END_REF][START_REF] Mazzini | Late Holocene palaeoenvironmental evolution of the Roman harbour of Portus, Italy[END_REF], Ostia-Rome [START_REF] Salomon | Geoarchaeology of the Roman port-city of Ostia: Fluvio-coastal mobility, urban development and resilience[END_REF], Neapolis-Naples [START_REF] Allevato | Pollen-wood analysis at the Neapolis harbour site (1 st -3 rd century AD, southern Italy) and its archaeobotanical implications[END_REF][START_REF] Donato | Development and decline of the ancient harbor of Neapolis[END_REF], Portus Pisanus-Pisa [START_REF] Kaniewski | Holocene evolution of Portus Pisanus, the lost harbour of Pisa[END_REF] or Pisa S. Rossore [START_REF] Mariotti Lippi | Pollen analysis of the ship site of Pisa San Rossore, Tuscany, Italy: the implications for catastrophic hydrological events and climatic change during the late Holocene[END_REF][START_REF] Bini | Palaeoenvironments and palaeotopography of a multilayered city during the Etruscan and Roman periods: early interaction of fluvial processes and urban growth at Pisa (Tuscany, Italy)[END_REF] have furnished key information on how human societies have shaped long-term landscape dynamics. While most studies have focused on the harbour phase, few analyses have focused on the post-harbour period, when structures were abandoned and anthropogenic pressures decreased or even disappeared (e.g. [START_REF] Marriner | The 'Ancient Harbour Parasequence': Anthropogenic forcing of the stratigraphic highstand record[END_REF][START_REF] Morhange | Eco-history of ancient Mediterranean harbours[END_REF].

Lying on the border between the Friuli plain and the edges of the Grado lagoon (Fig. 1), Aquileia, a Roman colony founded in 181 AD, and a UNESCO site since 1998, was located in a highly dynamic natural environment shaped by fluvial activity [START_REF] Arnaud-Fassetta | The site of Aquileia (northeastern Italy): example of fluvial geoarchaeology in a Mediterranean deltaic plain. Géomorphologie : relief, processus[END_REF].

Fresh, brackish and sea waters alternatively presented opportunities and risks for local communities and the site has been the object of long-term archaeological investigations and projects. During the Roman period, Aquileia was the focus of a communication network, consisting of terrestrial routes implemented by a system of natural and artificial waterways which provided access to the sea and maritime trade routes [START_REF] Cottica | Spunti per uno studio dell'interazione uomo e fiume in antico: il caso della sponda orientale del Natiso cum Turro ad Aquileia[END_REF]. The city soon became the most important trading centre (glass, wine, oil, olives, wool, gold and spices) of the northern Adriatic area [START_REF] Zaccaria | Moenibus at portu celeberrima Aquileia: Storia di una città[END_REF]. By the 1 st century AD, Aquileia had become a thriving centre, a major emporium in the Mediterranean and the Adriatic Sea, and a key site in commercial routes, lying at a crossroads for trade goods transiting from the Mediterranean Sea to the Danube provinces [START_REF] Carre | Les fouilles du port fluvial d'Aquilée[END_REF][START_REF] Maggi | Lo scavo di Canale Anfora[END_REF].

Furthermore, the site played a key role in Roman imperial political and military history: from the age of Marcus Aurelius onwards, it withstood several attacks until it fell under the siege of Attila in 452 AD, following the invasion of the Huns. Subsequently, Aquileia lost its importance for ~200 years before regaining its status as the region's principal city in the 7 th century AD with the rise of patriarchal rule. Aquileia gradually declined during the 15 th century due to the increasing power of the Republic of Venice [START_REF] Capulli | Ships of Aquileia. Underwater archaeological research on marine and inland routes of the Upper Adriatic Sea[END_REF].

Despite its important historical status, few environmental studies have been undertaken at Aquileia [START_REF] Arnaud-Fassetta | The site of Aquileia (northeastern Italy): example of fluvial geoarchaeology in a Mediterranean deltaic plain. Géomorphologie : relief, processus[END_REF], 2010;[START_REF] Siché | Cartographie hydrogéomorphologique et paléochenaux fluviatiles en milieux profondément modifiés par les sociétés. L'exemple du port fluvial antique d'Aquilée dans la plaine du Frioul (Italie septentrionale, Adriatique)[END_REF], one of the most significant fluvial harbours of the Roman Empire [START_REF] Carre | Les fouilles du port fluvial d'Aquilée[END_REF][START_REF] Cottica | Per uno studio integrato di uomo e ambiente ad Aquileia: alcune riflessioni a partire da recenti indagini archeologiche[END_REF][START_REF] Cottica | Spunti per uno studio dell'interazione uomo e fiume in antico: il caso della sponda orientale del Natiso cum Turro ad Aquileia[END_REF]. The geographical location of the city, however, makes it a key area to probe the environmental transitions that have occurred over the past centuries in the northern Adriatic coastal plain.

Here we investigate the post-Roman ecosystem dynamics of Aquileia to i) reconstruct how the environments have evolved during the last 1500 years, ii) investigate the major tipping points in ecosystem dynamics and iii) probe the environmental and climate forcing factors. Using a fossil record, we provide an in-depth ecological study of the area, reconstructing the major vegetation changes and their evolution through time. We also furnish a pollen-based climate reconstruction to identify the climate changes during the pre-industrial era [Late Antique Little Ice Age (LALIA), Medieval Climate Anomaly (MCA), and Pre-industrial Little Ice Age (LIA)].

Aquileia shows that harbours are key areas to understand how the Mediterranean environment was shaped by both human societies and natural factors during recent centuries.

Present-day vegetation

The area is located at the boundary of the Mediterranean belt. Though extensively transformed by human activities, nowadays, the vestiges of the fluvial harbour area are surrounded by patches of Quercus robur, Carpinus betulus, Crataegus monogyna, Ligustrum vulgare and Cornus sanguinea with man-made stands of Pinus pinea and Cupressus sempervirens along the old quays. Several species of Prunus (P. domestica, P. avium, P. mahaleb, P. spinosa) develop in the adjacent areas. The riparian vegetation is composed of rare patches of Salix alba, Populus alba and Alnus glutinosa along with a dense wet meadow composed of Filipendula ulmaria, Agrostis stolonifera, Phragmites sp. and a Magnocaricion (including Carex sp., Cyperus sp., Mentha aquatica and Sagittaria sp.).

Material and methods

The core

Biological indicators were extracted from the 450-cm Aquileia core III (AQ III) drilled in the ancient Roman harbour (45°46'29.068"N, 13°22'10.376"E, 1 m a.s.l.) during a field campaign in 2020 (Fig. 1). The grain-size analysis was performed with H2O2 treatment and removal of organic matter. Only the granular fraction was analysed. The lithology of the core AQ III is organized as follows: the bottom two meters (450-250 cm) comprise sandy deposits overlain by 142 cm of silt, organic silt and peat (250-108 cm). The upper levels correspond to a siltysandy matrix (108-0 cm; Fig. 2). No evidence of stratigraphic hiatuses was observed in the core or laboratory data. The core was sampled every 2 cm on average. Bioindicators are only wellpreserved in the silt-organic-peat layers.

Chronology

The chronology of the core AQ III (Fig. 2) is based on seven accelerator mass spectrometry (AMS) 14 C dates performed on short-lived samples (small leaves) at 110 cm (200±30 14 C yr BP), 126 cm (270±30 14 C yr BP), 156 cm (380±30 14 C yr BP), 166 cm (470±30 14 C yr BP), 204 cm (1260±30 14 C yr BP), 232 cm (1490±30 14 C yr BP), and 246 cm(1540±30 14 C yr BP).

Aquatic remains were excluded, and all of the botanical macro-remains sent to Beta Analytic (Florida) comprised small leaves from deciduous vegetation. The 14 C dates were calibrated to 2 sigma (σ) using Calib Rev 8.0.1 (IntCal20). All the calibrated ages are denoted in AD, consistent with the historical data. The age-depth model is based on a linear interpolation between each of the intercepts (and the 2σ). The interpolation is compared with a linear model and a polynomial model in order to estimate the difference between the two signals (Fig. 2).

Sediment compaction during coring was calculated in order to obtain the real depth of each deposit. Subsequently, the sedimentation rate between each 14 C date was calculated (Fig. 2).

Because samples were taken at regular intervals along the 250-cm sediment column, the chronological resolution is directly dependent on the sedimentation rate (Fig. 2).

Biological data

Samples from the core AQ III were prepared for pollen analysis using standard procedures for clay-silt samples (fully detailed in [START_REF] Faegri | Textbook of pollen analysis[END_REF]. Pollen grains were counted under x400 and x1000 magnification using an Olympus microscope. Pollen frequencies (expressed as percentages) are based on the terrestrial pollen sum, excluding local helophytes, macrophytes and spores of non-vascular cryptogams. Aquatic taxa frequencies were calculated by adding the local helophytes-macrophytes to the terrestrial pollen sum. The mean pollen sum was 507±145 pollen grains, with a minimum of 240 pollen grains (there were only five samples with less than 300 pollen grains). The median value was 547 pollen grains, with a 25 th percentile of 383 pollen grains and a 75 th percentile of 627 pollen grains. The mean number of taxa was 59±5 taxa, with a minimum of 42 taxa. The median value was 60 taxa, with a 25 th percentile of 58 taxa, and a 75 th percentile of 63 taxa. The palynological data are presented as a detailed pollen diagram in Fig. 3. The river floodplain vegetation is shown as a second detailed pollen diagram (Fig. 4).

To support the identification of the main environmental changes along the core, sixteen samples were also collected for qualitative analysis of the ostracod fauna with a resolution of <15 cm.

Samples were prepared adopting the standard methodology (e.g. [START_REF] Amorosi | Middle to late Holocene environmental evolution of the Pisa coastal plain (Tuscany, Italy) and early human settlement[END_REF] and the taxonomic identification was based on [START_REF] Henderson | Freshwater ostracods[END_REF].

Statistical analyses

Statistical analyses were performed using Xl-Stat 2019 (https://www.xlstat.com), R 4.1.0 (R Core Team, 2021) and the software package PAST 4.08 [START_REF] Hammer | PAST: paleontological statistics software package for education and data analysis[END_REF]. All the long-term trends were calculated using a polynomial model (order 5 -lowest standard deviation, Pvalue < 0.001; the R 2 is mentioned on each graph). Pollen data were studied using a cluster analysis (descending type; Fig. 5) to calculate a dendrogram, using branches as ecological distances between groups of taxa. The test was performed using Paired group as algorithm and Correlation as the similarity measure. Each cluster was summed to generate pollen-derived vegetation patterns (PdVs) and assigned to a potential location, from the dryland to the floodplain, referring to modern patches of vegetation (Fig. 5). The riparian forest (fen trees) was added to the marshland plants (macrophytes and helophytes) to create a curve reflecting the floodplain dynamics (Fig. 6A). A principal component analysis (PCA) was subsequently performed to test the ordination of samples by assessing major changes in the PdV-scores. The first axis (PCA-Axis 1), which carries the maximum variance (77% of total variance), was extracted. The PdVs influenced by the floodplain are loaded by the positive PCA-Axis 1 scores while the PdVs located beyond the influence of the river are loaded by the negative scores. A regular chronological interpolation (10-yr) was applied to the dataset in order to reduce biases linked to chronological gaps, particularly at the bootstrap level. We then transformed the PCAAxis1 scores using a Loess smoothing (with a LOWESS algorithm) and performed a bootstrap to estimate a 95% confidence band based on 1000 random replicates. The Loess curve and the 95% confidence band were used as a proxy for the floodplain dynamics (Fig. 6B). The cultivated species and pasture activities were added to test their occurrences in a context of long-term floodplain changes (Fig. 6C).

Variations in the floodplain (fen trees, helophytes and macrophytes) were further contrasted with the dynamic of the dry floodplain scrubs, using both a synthetic diagram and a Kernel density 2D model (with Gaussian as function; Fig. 7A). The two time-series were then tested using a cross-correlation (Pvalue < 0.001; the R 2 value is indicated on Fig. 7A). The floodplain dynamics (fen trees and helophytes) were subsequently compared to the dryland components (pine-oak woodland cluster), using both a synthetic diagram and a Generalized Linear Model (Fig. 7B). The two time-series were also tested using a cross-correlation (Pvalue < 0.001; the R 2 value is indicated on Fig. 7B). A Gradient Species Packing test was finally added to identify the distribution of each cluster, using the floodplain as a stationary reference (Fig. 7C).

Agriculture, pasture activities and fire

Agriculture activities are mainly based on cereals (Poaceae cerealia) with a diversified arboriculture comprising grapevines (Vitis vinifera), olive trees (Olea europaea), walnut trees (Juglans regia) and diverse Rosaceae (Prunus sp.). Pasture activities are identified by various genera such as knapweed (Centaurea), plantain (Plantago), knotweed (Polygonum) and sorrel (Rumex). To identify human impacts, a pasture activities/agriculture ratio was calculated (Fig. 8A). All the long-term trends were calculated using a polynomial model (order 5, Pvalue < 0.001).

The fire history of Aquileia was reconstructed by counting charcoal particles (50-200 μm) on the pollen slides. The 50 μm size criterion was chosen to avoid confusion between microscopic charcoal fragments and opaque minerals, which are typically < 50 μm [START_REF] Parshall | Fire on the New England landscape: regional and temporal variation, cultural and environmental controls[END_REF][START_REF] Pederson | Medieval Warming, Little Ice Age, and European impact on the environment during the last millennium in the lower Hudson Valley[END_REF]. We adopted the 200-300 items (sum of charcoal particles) defined by [START_REF] Finsinger | Minimum count sums for charcoal concentration estimates in pollen slides: accuracy and potential errors[END_REF]. A first charcoal concentration curve (fragments per cm -3 ) was drawn according to a linear age-scale (Fig. 8A). The fire activity was compared with the ratio pasture activities/agriculture, with the dynamics of the mixed oak forest (Fig. 8B), and with agricultural activities (Fig. 8C). Charcoal concentrations were converted into charcoal accumulation rates (CHAR) based on the sedimentation rate estimated by the age-depth model (Fig. 8D). These analyses were performed in R using R-PaleoAnomalies (https://github.com/wfinsinger/R-PaleoAnomalies), a code which builds on CharAnalysis. The analysis comprised the following steps: i) The CHAR record was resampled to equal time intervals (20-yr) using a window corresponding to the mean sedimentation rate of the record, ii) the data were log-transformed to homogenize variance, iii) the "Background Component" was estimated with a local polynomial (100-yr moving window), facilitating the calculation of the "Peak Component", and iv) fire episodes were detected using a threshold within the range of values with the lowest sensitivity to the number of peaks detected, derived by plotting a frequency distribution histogram of the peak component. The process is fully described in the literature [START_REF] Long | A 9000-year fire history from the Oregon Coast Range, based on a high-resolution charcoal study[END_REF][START_REF] Higuera | Frequent fires in ancient shrub tundra: implications of paleorecords for arctic environmental change[END_REF][START_REF] Higuera | Vegetation mediated the impacts of postglacial climate change on fire regimes in the southcentral Brooks Range, Alaska[END_REF][START_REF] Doyen | Climate and biomass control on fire activity during the late-glacial/early-Holocene transition in temperate ecosystems of the upper Rhone valley (France)[END_REF].

Pollen-derived climate reconstruction

The pollen-based model used to reconstruct climate variables (Fig. 9) has been described extensively in the literature [START_REF] Cheddadi | Climate change since the last glacial period in Lebanon and the persistence of Mediterranean species[END_REF]Cheddadi et al., 2016[START_REF] Cheddadi | Microrefugia, climate change, and conservation of Cedrus atlantica in the Rif Mountains, Morocco[END_REF] and has been shown to provide coherent climate reconstructions (e.g., [START_REF] Kaniewski | Cold and dry outbreaks in the eastern Mediterranean 3200 years ago[END_REF][START_REF] Kaniewski | Recent anthropogenic climate change exceeds the rate and magnitude of natural Holocene variability on the Balearic Islands[END_REF][START_REF] Kaniewski | Coastal submersions in the north-eastern Adriatic during the last 5200 years[END_REF]. Past climate variables were derived from the Aquilea fossil pollen record. We assigned the pollen taxa of Aquilea to modern plant species based on our knowledge in palynology, botany, and ecology. Edaphic, aquatic, human-related, or less than 1% of the total pollen taxa were not considered in the climate reconstruction. For each fossil sample, we used the weighted median of each pollen taxon identified in individual fossil samples, using the percentages of pollen as weights to calculate the median of all pollen taxa. The uncertainty of the reconstruction for each fossil sample is estimated by omitting one taxon from each iterative estimate, as many times as the number of taxa in the fossil sample. This quantification method assumes that the modern climatic range of a species encompasses its range during the studied time span and that the closer the best climatic value is to the area where it is most abundant, the more abundant it is [START_REF] Cheddadi | Early Holocene Greening of the Sahara requires Mediterranean winter rainfall[END_REF]. Helophytes, macrophytes and spores of nonvascular cryptogams were excluded from the matrix. All reconstructions are shown with their estimated errors (Fig. 9), and with their long-term trends based on a polynomial model (order 5, Pvalue < 0.001; the R 2 is shown on each graph). The temperature anomalies are based on the annual average for the control period 1961-1990 (13.04°C; Fig. 9A). Temperature anomalies were first compared with the loess-smoothed PCA-Axis 1 (Fig. 9B) and the spring-summer precipitation (Fig. 9C). The 1850-2019 CRUTEM data [START_REF] Osborn | The CRUTEM4 land-surface air temperature data set: construction, previous versions and dissemination via Google Earth[END_REF] were then added to test the reliability of the temperature reconstruction from core AQ III (Fig. 9D). A sinusoidal model was applied (Pvalue < 0.001) to link data from core AQ III with CRUTEM4 (Fig. 9D).

Results

Lithology and chronology

The lower part of the core AQ III, corresponding to an abandoned channel of Roman age, was completely void of bioindicators. The sandy matrix was not conducive to the conservation of ostracods, pollen and other non-palynomorphs. The overlying levels, which correspond to silt, organic silt and peat deposits (Fig. 2), were found to be rich in palynomorphs, with an average of 507±145 pollen grains counted per sample and an average concentration of 2027±585 pollen per cm -3 . The median value is 2186 pollen per cm -3 with a 25 th percentile of 1531 pollen per cm -3 and a 75 th percentile of 2506 pollen per cm -3 . Within this interval, hypohaline ostracods are locally found although represented by few valves. The sedimentation rate (SR) varies according to the lithology. While the lower silty layers have a SR of 1.7 to 2 mm per year, the organic silt and peat deposits are characterized by lower rates of 0.5 mm per year (Fig. 2). The upper silty layers are characterized by a SR of 1.25 to 2 mm per year, close to the values recorded for the lower layers. The uppermost level corresponds to a silty-sandy matrix completely sterile of bioindicators. The material used for 14 C dating was only found in the siltorganic and silt-peat layers. No macro-remains were observed in the lower sandy matrix and upper silty-sand deposits. The core AQ III is framed by a chronology spanning from 510±85

AD to 1770±40 AD (Fig. 2), or from the Ostrogothic Kingdom to the Austrian Empire.

Ecosystem dynamics

Since 500 AD, the evolution of vegetation surrounding the waterlogged area developed after the abandonment of the Roman fluvial harbour is detailed in two pollen diagrams with the first focused on terrestrial plants (Fig. 3) and the second on floodplain vegetation (Fig. 4). The pollen dataset was categorized into 9 PdVs grouped into 4 assemblages (Dryland, Dry floodplain, Wetland, Floodplain) using a cluster analysis (Fig. 5). The PdVs were then converted into markers of hydrological shifts using the floodplain vegetation (Fig. 6A) and a Loess PCA-Axis 1 (Fig. 6B) as proxies. The outcomes suggest a humid environment until 660±55 AD, with a well-developed floodplain (macrophytes and helophytes) surrounded by fen trees (Alnus, Betula, Populus and Salix) and by a wet meadow (dominated by Poaceae, Ranunculaceae and Campanulaceae). Most of the macrophytes and helophytes suggest the presence of stagnant waters or a river branch with a very low streamflow, creating a swamp in the city. This interpretation is consistent with the occurrence of ostracods typical of shallow hypohaline environments with stagnant waters (Pseudocandona and Cyclocypris species). Two spikes at 565-590±55 AD may be consistent with floodplain high-stands (Fig. 6A). Agriculture was weakly developed during this first phase (Fig. 6C). After 660±55 AD, the floodplain regressed significantly, reaching its first low level at 950±45 AD and a second at 1200-1230±35 AD (Fig. 6). The macrophytes and helophytes were gradually replaced by dry floodplain scrubs (Fig. 7A)

and by a pine-oak woodland (Fig. 7B), suggesting the development of drier environmental conditions. Accordingly, the disappearance of ostracods points to a marked decrease in the water table. A reliable linear relationship (R 2 =0.66; Pvalue < 0.001) was established between the loss of fen trees and helophytes, and the increase of pine-oak woodland, implying an inverse relationship. During this period (660±55 to 1230±35 AD), agriculture and pastoral activities increased significantly (Fig. 6C). From 1230±35 to 1770±40 AD, the area oscillates between dry and wet conditions, suggesting unstable environmental conditions. The drier phases are supported by xeric components (e.g. Chenopodiaceae, Artemisia, Asteraceae, Quercus evergreen). The wet phases, consistent with a rise in the floodplain's water table, are attested by an increase in helophytes and wet meadow plants (Fig. 4), and the local reappearance of a hypohaline ostracod fauna. After 660±55 AD, and the drying of the swamp, wet phases were recorded at 1000-1050±40, 1250-1310±35, 1380-1420±25, 1480-1500±40, and 1590-1650±30 AD (Fig. 6A-B). After 1200±35 AD, agricultural and pastoral activities declined sharply, then fluctuated, before regaining in importance after 1710±35 AD (Fig. 6C).

Agricultural practices, pastoral activities and fire

The agricultural practices (agri-and arboriculture) and pastoral activities were first depicted as long-term dynamics (Fig. 6C) and then shown as ratios (Fig. 8A). During the last 1500 years, agriculture has never exceeded 15% of the pollen sum (with a maximum of 6% for the cereals)

and the pastoral plants never rise above 12%. Both were probably more developed during the Roman period and the heydays of the city. Agriculture is mainly recorded from 620±55 to 1230±35 AD, during the patriarchal rule of the city, and then after 1710±35 AD, during the Austrian Empire (Fig. 8A). Fire activity meshes inversely with the mixed oak forest phase (Fig. 8A). The mixed oak forest regressed most significantly when fire return intervals were high, and inversely developed when slash-and-burn agriculture declined (Fig. 8B). When both agricultural and pastoral practices are combined, there is a clear link between anthropogenic activities and ignitions, suggesting human-induced fires in the area (Fig. 8C and8D). From 1230±35 to 1430±35 AD, a period marked by general decline in anthropogenic and cultivated taxa, the pastoral activities/agriculture ratio shows a prevalence of pastoral activities over agriculture (Fig. 8A). After 1710±35 AD, a declining ratio is consistent with a renewed increase in agricultural activities (Fig. 8A), probably slash-and-burn agriculture as fire indicators also increase (Fig. 8D).

Climate reconstruction

The reconstruction of temperature anomalies since 500 AD depicts colder periods centred on the LALIA, the Oort minimum, and the different episodes of the LIA (Wolf minimum, Spörer minimum, Maunder minimum; Fig. 9A; [START_REF] Marriner | Anthropocene tipping point reverses long-term Holocene cooling of the Mediterranean Sea: A meta-analysis of the basin's Sea Surface Temperature records[END_REF]. The colder phase corresponds to the LALIA, locally recorded from 535±75 to 660±55 AD, with a minimum of -3.5±0.17°C at 590±55 AD and an average of -2.04±0.17°C for the whole phase. The second cooler period corresponds to the Maunder episode, with a minimum of -1.4±0.14°C and an average of 0.78±0.15°C for the whole phase. The cold periods are correlated with greater humidity and the expansion of the swamp area, as highlighted by the Loess PCA-Axis 1. An exception is the Maunder minimum, during which the area remained dry (Fig. 9B). Focusing on spring-summer precipitation, the drier phases are centred on the Medieval Climate Anomaly (average 291±2 mm), and on the period lasting from the end of the Wolf minimum to the onset of the Maunder minimum (average 287±2 mm; Fig. 9C). While precipitation increases are concomitant with the development of the swamp area, the Maunder minimum is characterized by higher rainfall but a dry context (Fig. 9B).

When the reconstructed temperature anomalies are compared and contrasted with the CRUTEM4 dataset (1850-2019 AD; [START_REF] Jones | Hemispheric and large-scale land-surface air temperature variations: An extensive revision and an update to 2010[END_REF][START_REF] Osborn | The CRUTEM4 land-surface air temperature data set: construction, previous versions and dissemination via Google Earth[END_REF], there is a clear link between the record from Aquileia and the climatic evolution of the northern Adriatic. The anomalies of -1.09±0.14°C at 1770±40 AD (last point recorded in the core AQ III) and of 1.13°C at 1850 AD (first data available in the CRUTEM4 dataset) are consistent. The sinusoidal model (Pvalue < 0.001) shows no discontinuity between the two curves (Fig. 9D), suggesting that core AQ III and CRUTEM4 are comparable. It also appears that the MCA was not warmer than the 20 th century AD in the northern Adriatic, but cooler compared to the 21 st century AD.

Discussion

The fluvial harbour of Aquileia offers a rare record of the environmental shifts that have occurred on the northern Adriatic coastal plain since 500 AD. In the Gulf of Venice, the late Holocene ecosystem dynamics have been poorly studied, especially the last two millennia.

Previous work has focused on the dynamics of salt-marsh plant communities [START_REF] Miola | Holocene salt marsh plant communities in the North Adriatic coastal plain (Italy) as reflected by pollen, non-pollen palynomorphs and plant macrofossil analyses[END_REF], on long-term sequences covering the Pleistocene-Holocene (e.g. [START_REF] Bortolami | Land, sea and climate in the northern Adriatic region during late Pleistocene and Holocene[END_REF][START_REF] Bondesan | Paleoenvironmental reconstruction from LGM to historical time in the lower coastal plain of the Piave river. Preliminary pollen analysis on a 20 m core of lagoon and fluvial sediments[END_REF][START_REF] Masari | The environment of Venice area in the past two million years[END_REF][START_REF] Pini | Pollen stratigraphy, vegetation and climate history of the last 215 ka in the Azzano Decimo core (plain of Friuli, north-eastern Italy)[END_REF][START_REF] Donnici | Evidence of climatic changes in the Venetian Coastal Plain (northern Italy) during the last 40,000 years[END_REF], on grasslands covering the eastern part of the northern Adriatic [START_REF] Kaligarič | Recent vegetation history of the north Adriatic grasslands: expansion and decay of an anthropogenic habitat[END_REF], and on sequences from the Istrian Peninsula, Croatia (e.g. [START_REF] Kaniewski | Solar pacing of storm surges, coastal flooding and agricultural losses in the Central Mediterranean[END_REF][START_REF] Kaniewski | Coastal submersions in the north-eastern Adriatic during the last 5200 years[END_REF].

A wet and cold LA Little Ice Age

From 508±80 to 660±55 AD, Aquileia was wet, marked by significant occurrences of fen trees, helophytes, macrophytes and by an expanded wet meadow. The floodplain that hosted the city of Aquileia was well-developed during the Ostrogothic Kingdom and the early phase of the Kingdom of the Lombards. The zone was a swamp area, characterized by low-energy fluvial dynamic in the center surrounded by stagnant water. The background vegetation is mainly composed of a mixed oak forest and by short incursions of dry floodplain scrubs. The upper part of the marine core RF 93-30, from the western flank of the Adriatic [START_REF] Oldfield | A high resolution late Holocene palaeo environmental record from the central Adriatic Sea[END_REF] depicts a similar vegetation pattern, with a dominant mixed oak forest and a strong presence of helophytes (Cyperaceae). At Aquileia, two major wet peaks at 565-590±55 AD may suggest a floodplain high stand (Fig. 6A). Similar conclusions were reached by a geomorphological study of the Natiso River [START_REF] Siché | Migrations et métamorphoses historiques des fleuves torrentiels sur leur delta et leurs impacts sur les implantations portuaires antiques. L'exemple de l'hydrosystème Torre Natisone Isonzo sur la plaine côtière d'Aquilée (Méditerranée nordoccidentale)[END_REF]. As agriculture and pastural activities were weakly developed in the northern Adriatic (Fig. 6C), the environmental variations seem to be mainly related to climate. The LALIA is defined as a cold episode which lasted from 536 to 660 AD. Some studies suggest that it exceeded the LIA in severity [START_REF] Büntgen | Cooling and societal change during the Late Antique Little Ice Age from 536 to around 660 AD[END_REF], with a major impact in Italy [START_REF] Neukom | No evidence for globally coherent warm and cold periods over the preindustrial Common Era[END_REF][START_REF] Peregrine | Climate and social change at the start of the Late Antique Little Ice Age[END_REF]. The reconstructed temperatures at Aquileia show a similar pattern, with an average of -2.04±0.17°C for the LALIA, and an average of 0.78±0.15°C for the coldest event of the LIA (Maunder minimum). The minima follow the same trend, with a value of -3.5±0.17°C at 590±55 AD (LALIA) and -1.4±0.14°C at 1642±30 AD (LIA; Fig. 9A). The LALIA at Aquileia was, on average, approximately -1.26°C colder than the LIA. A wetter and colder climate, as shown by the core AQ III (Fig. 9C), has also been elucidated in the Gulf of Gaeta [START_REF] Di Rita | Late Holocene forest dynamics in the Gulf of Gaeta (central Mediterranean) in relation to NAO variability and human impact[END_REF], in Sicily [START_REF] Sadori | Climate, environment and society in southern Italy during the last 2000 years. A review of the environmental, historical and archaeological evidence[END_REF] and in coastal Croatia where spring-summer precipitation increased significantly during this phase [START_REF] Kaniewski | Coastal submersions in the north-eastern Adriatic during the last 5200 years[END_REF]. At other sites, this period is inversely associated with a marked human impact on the environment (e.g. Lago Lungo; [START_REF] Mensing | 2700 years of Mediterranean environmental change in central Italy: a synthesis of sedimentary and cultural records to interpret past impacts of climate on society[END_REF]core ND2;[START_REF] Michelangeli | Three millennia of vegetation, land-use and climate change in SE Sicily[END_REF]. The core AQ III suggests a cold and wet northern Adriatic from 508±80 to 660±55 AD consistent with Mediterranean Sea Surface Temperatures at this time [START_REF] Marriner | Anthropocene tipping point reverses long-term Holocene cooling of the Mediterranean Sea: A meta-analysis of the basin's Sea Surface Temperature records[END_REF].

A shift towards warmer and drier conditions

From 660±55 to 1230±55 AD, the area changed radically from a floodplain to a drier environment as shown by the Loess PCA-Axis1 (Fig. 6B), with an important extension of the scrubs on desiccated banks, and of the pine-oak woodland in adjacent areas. The Natiso River was transformed into a minor waterway with reduced flow. One of the main changes since 660±55 AD was the development of agriculture during the patriarchal rule of the city. Fire has long been an element of pressure on ecosystems and an agent associated with agricultural activities in Italy (e.g. [START_REF] Sadori | Charcoal analysis, a method to study vegetation and climate of the Holocene: The case of Lago di Pergusa (Sicily, Italy)[END_REF][START_REF] Vannière | Climate versus human-driven fire regimes in Mediterranean landscapes: the Holocene record of Lago dell'Accesa (Tuscany, Italy)[END_REF]. At Aquileia, the same pattern was found, but agriculture was however limited as taxa associated with anthropogenic activities never exceed 15% of the pollen sum during this period. The olive trees never reached values equivalent to those recorded at the Lago Alimini Piccolo in southern Italy [START_REF] Di Rita | Holocene drought, deforestation and evergreen vegetation development in the central Mediterranean: a 5500-year record from Lago Alimini Piccolo, Apulia, southeast Italy[END_REF] and are closer to the values identified at Lago Lungo [START_REF] Mensing | 2700 years of Mediterranean environmental change in central Italy: a synthesis of sedimentary and cultural records to interpret past impacts of climate on society[END_REF] and in the core ND2 [START_REF] Michelangeli | Three millennia of vegetation, land-use and climate change in SE Sicily[END_REF]. Similarly, fire events are mainly recorded during the reintroduction of cultivated species (Fig. 8D), a short period extending from 660±55 to 750±55 AD and after 1700 AD. The development of slash-and-burn agriculture, while limited, seems to have mainly affected the mixed oak forest, freeing up ecological niches colonized by the pine-oak woodland (Fig. 7B). As agriculture was not dominant in the area, the variations recorded by the core AQ III derive from a mixed climate-anthropogenic signal, with an environmental forcing more influential than anthropogenic factors as suggested by the occurrence of the main climatic phases of the pre-industrial era. Focusing on the climate parameters, it appears that the area was warmer after the LALIA, with a first cooler period corresponding to the Oort minimum (from 1000±40 to 1050±40 AD; Fig. 9A). The Oort minimum, defined as a period of low solar activity [START_REF] Yiou | Stability of weather regimes during the last millennium from climate simulations[END_REF], does not seem to have been a major event in the northern Adriatic, with an average of -0.2±0.14°C recorded at Aquileia. The weak influence of the Oort minimum in the Gulf of Venice has already been highlighted with regards to sea-storm activity [START_REF] Camuffo | Sea storms in the Adriatic Sea and the western Mediterranean during the last millennium[END_REF]. After 1050±40 AD, the area around Aquileia became warmer, with temperature anomalies peaking at 0.4±0.15°C at 1135±35 AD. While the MCA (core period: 1000-1200 AD) is characterized by anomalously warm conditions in some (but not all) regions [START_REF] Mann | Global signatures and dynamical origins of the Little Ice Age and Medieval Climate Anomaly[END_REF][START_REF] Lüning | The Medieval Climate Anomaly in the Mediterranean Region[END_REF], in the northern Adriatic, the average temperature anomaly for the whole episode is 0.08±0.15°C, close to the control period 1961-1990 AD (13.04°C) AD (13.04°C). This observation is corroborated by a European study which revealed that temperatures during the MCA were probably as warm as the 20 th century AD [START_REF] Büntgen | European tree-ring data and the Medieval Climate Anomaly[END_REF]. The MCA at Aquileia ended at 1230-1250±35 AD. AQ III shows a warm and dry northern Adriatic from 660±55 to 1230-1250±35 AD, with a low impact of the Oort minimum and a local, limited, increase of anthropogenic pressures on the environment.

A succession of cold and warm phases

From 1250±35 to 1770±40 AD, the swamp area oscillated between wet and dry phases. The dry phases are attested by the occurrence of xeric components within and beyond the floodplain while the wet episodes are mainly attested by an increase in helophytes and wet meadow plants.

This suggests that the flow of Natiso River has oscillated markedly over time (Fig. 6A).

Although the wet phases are clearly identified by the Loess PCA-Axis1, the waterbody never reached the level recorded during the LALIA. Following the end of patriarchal rule (in 1236 AD), agricultural, pastoral and fire activities decreased markedly. All the anthropogenic pressures on ecosystems seem to have been reduced until 1710±35 AD and renewed after this period (Fig. 6C). The signal is thus mainly modulated by environmental factors, namely climate variations. The first cold and wet phase occurred from 1250±35 to 1335±30 AD and meshes with the Wolf minimum, a Grand Solar Minimum spanning the period 1270-1350 AD [START_REF] Steinhilber | Total solar irradiance during the Holocene[END_REF][START_REF] Usoskin | Solar activity during the Holocene: the Hallstatt Cycle and its consequence for Grand Minima and Maxima[END_REF]Fogtmann-Schulz et al., 2021). At Aquileia, the temperature anomaly recorded during this phase is on average -0.32±0.16°C. The area became colder but also wetter, with a clear positive deviation of the spring-summer precipitations, similar to the trends observed in coastal Croatia [START_REF] Kaniewski | Coastal submersions in the north-eastern Adriatic during the last 5200 years[END_REF]. The second cold and wet phase occurred from 1460±35 to 1500±40 AD and corresponds to the central point of the Spörer minimum that spans the period 1390-1550 AD [START_REF] Fogtmann-Schulz | Variations in Solar Activity Across the Spörer Minimum based on Radiocarbon in Danish Oak[END_REF][START_REF] Marriner | Anthropocene tipping point reverses long-term Holocene cooling of the Mediterranean Sea: A meta-analysis of the basin's Sea Surface Temperature records[END_REF]. This episode is weakly recorded at Aquileia, with a maximum anomaly of 0.09±0.14°C. Inversely, this phase is characterized by an important positive deviation in springsummer precipitation, correlated with a catastrophic flood event in the nearby Isonzo River at 1480 AD [START_REF] Siché | Migrations et métamorphoses historiques des fleuves torrentiels sur leur delta et leurs impacts sur les implantations portuaires antiques. L'exemple de l'hydrosystème Torre Natisone Isonzo sur la plaine côtière d'Aquilée (Méditerranée nordoccidentale)[END_REF] and a significant increase in precipitation in coastal Croatia [START_REF] Kaniewski | Coastal submersions in the north-eastern Adriatic during the last 5200 years[END_REF]. This event is also correlated with a maximum frequency of sea storms in the Gulf of Venice [START_REF] Camuffo | Sea storms in the Adriatic Sea and the western Mediterranean during the last millennium[END_REF]. The last cold and wet episode recorded in AQ III fits with the Maunder minimum (1645-1715 AD; Shindell et al., 2001;[START_REF] Usoskin | The Maunder minimum (1645-1715) was indeed a grand minimum: A reassessment of multiple datasets[END_REF] that locally spans the period 1620±30 to 1690±30 AD. On average, the temperature anomaly is 0.78±0.15°C with a minimum of -1.4±0.14°C (Fig. 9A). The Maunder minimum is the second coldest episode and also the second wettest event, after the LALIA. The warmer phase is recorded from 1520±40 to 1600±30 AD, with an average annual temperature anomaly of 0.71±0.15°C. This warm period is correlated with a major increase in total solar irradiance during a phase of low solar activity [START_REF] Steinhilber | Total solar irradiance during the Holocene[END_REF], and higher winter and summer temperatures in coastal Croatia [START_REF] Kaniewski | Coastal submersions in the north-eastern Adriatic during the last 5200 years[END_REF]. The Eurasian summer temperature variability reconstructed from tree rings also depicts anomalies close to the MCA during this period [START_REF] Büntgen | Cooling and societal change during the Late Antique Little Ice Age from 536 to around 660 AD[END_REF]. When the anthropogenic activities resumed at 1710±35 AD under the Austrian Empire (Fig. 8C), the climate conditions were warm and still wet, facilitating agriculture and arboriculture. AQ III points to unstable environmental conditions in the northern Adriatic from 1250±35 to 1770±40 AD, with a maximum impact during the Maunder minimum.

Even though AQ III is chronologically correlated with solar irradiance, a direct causal relationship is not suggested here because the cold LIA periods began earlier than the decline in solar activity and the Grand Solar Minimum.

Conclusions

The port city of Aquileia has provided a rare reconstruction of climatic and anthropogenic pressures on ecosystems from the end of the Roman Empire to the Austrian Empire in the northern Adriatic. The abandoned fluvial harbour evolved into a swamp, recording the major tipping points in the ecosystem dynamics, mainly due to climate pressures and, secondly, to anthropogenic activities. The main phases recorded by the ecosystems correspond to the LALIA and the Maunder minimum, with the development of cold and wet conditions. Inversely, the positive temperature anomalies reconstructed for the MCA suggest that this episode is close to temperatures recorded for the 20 th century AD but, however, they differ from those of the 21 st century AD. Anthropogenic activities, while diverse, seem to have been a secondary pressure at Aquileia. Their development is chronologically limited (late 7 th to early 13 th centuries AD)

and their amplitude not strong enough to supersede the other environmental pressures. The long-term ecosystem changes that may correspond to human pressures on the rural environment are probably masked by the climate dynamics. Aquileia suggests that the environmental pressures have been key in shaping the landscapes of the climate-sensitive northern Adriatic.

Figures

Figure 1 .

 1 Figure 1. Geographical location of the study area in Italy. The ancient harbour of Aquileia is marked in red. The city of Venice is indicated on the map by a green circle.

Figure 2 .

 2 Figure 2. Focus on the core AQ III and the radiocarbon chronology. The lithology of the core is detailed according to depth. The sedimentation rates are shown in mm per year. The radiocarbon dates are depicted as intercepts and 2-sigma calibrations (95% of probability). The age model is compared and contrasted with linear (purple line) and polynomial (orange line) regressions.

Figure 3 .

 3 Figure 3. Detailed pollen diagram (woody and herbaceous vegetation) for the period 5001800 AD. The taxa are plotted on a linear age-scale. The vertical lines indicate the main transitions.

Figure 4 .

 4 Figure 4. Detailed pollen diagram of the Aquileia river floodplain for the period 500-1800 AD. The taxa are grouped into PdVs and plotted on a linear age-scale. The vertical lines indicate the main transitions.

Figure 5 .

 5 Figure 5. Pollen-based ecological clusters from Aquileia for the period 500-1800 AD. A cluster analysis (Paired group as algorithm, Correlation as similarity measure) was used to define the ecological assemblages. Each cluster was summed to create the PdVs. The PCA loadings of each PdV are indicated on the graph.

Figure 6 .

 6 Figure 6. Environmental dynamics of the floodplain. (A) Evolution of the floodplain depicted by the riparian forest and marshland. The occurrence of hypohaline ostracods is also shown. (B) Loess PCA-Axis 1 (77% of total variance) shown with the 95% confidence band. (C) Pastoral activities alone (dark green) and pastoral activities with agriculture (light green) are plotted on a linear age-scale. All the long-term trends are highlighted by polynomial models (order 5, Pvalue < 0.001).

Figure 7 .

 7 Figure 7. Floodplain dynamics compared and contrasted with dry components. (A) Variations in the floodplain contrasted with the dynamics of the dry floodplain scrubs based on a synthetic diagram and a 2D Kernel-density model. The R 2 of the cross-correlation is indicated on the Kernel graph (B) Floodplain dynamic compared to the dryland components based on a synthetic diagram and a Generalized Linear Model. The R 2 of the cross-correlation is indicated on the synthetic diagram (C) Gradient Species packing test showing the distribution of each cluster taking the floodplain as a stationary reference. All the long-term trends are highlighted by polynomial models (order 5, Pvalue < 0.001).

Figure 8 .

 8 Figure 8. Agriculture, pasture and fire activity. (A) Ratio pasture activities/agriculture and fire dynamics plotted on a linear age-scale. (B) Mixed oak forest dynamics compared and contrasted with fire activity (shown as a shaded curve). (C) Agro-pastoral activities contrasted with the fire activity (shown as a shaded curve). (D) CHAR analysis plotted on a linear agescale.The red circles with white crosses denote the fire episodes. All the long-term trends are depicted by polynomial models (order 5, Pvalue < 0.001).

Figure 9 .

 9 Figure 9. Pollen-based climate reconstruction. (A) Temperature anomalies (with the standard deviations) plotted on a linear age-scale. (B) Loess PCA-Axis 1 (77% of the total variance) shown with the 95% confidence band. (C) Spring-summer precipitation (with the standard deviations) plotted on a linear age-scale. All the long-term trends are depicted by polynomial models (order 5, Pvalue < 0.001). (D) Temperature anomalies compared and contrasted with the 1850-2019 CRUTEM4 data[START_REF] Osborn | The CRUTEM4 land-surface air temperature data set: construction, previous versions and dissemination via Google Earth[END_REF]. The sinusoidal model (Pvalue < 0.001) is depicted as a purple line.

  Figure 1

Figure 7 Figure 8 Figure 9

 789 Figure 7

Acknowledgments

We wish to thank Daniele Pittaro for his support during the fieldwork campaign at Aquileia. The paper presents results of post-excavation laboratory analyses carried out as part of a wider project of archaeological and palaeoenvironmental research undertaken in the fluvial harbour of Aquileia (ex Fondo Sandrigo and adjacent areas) by the University of Venice Ca' Foscari, headed by D. Cottica, and under the auspices of the Italian Ministry of Culture in collaboration with the Soprintendenza Archeologia, belle arti e paesaggio del Friuli Venezia Giulia (Sabap FVG).

Data availability

All data needed to evaluate the conclusions in the paper are present in the paper and/or the

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Supplementary Materials (Raw data file).

Funding

Support was provided by the MITI CNRS "Evénements rares", AQUASANMARCO program, and the EPHE-DRI program. On-site research was funded by the Ca' Foscari University of Venice in collaboration with the Fondazione Aquileia.