Matthieu Furet
email: matthieu.furet@ls2n.fr

Philippe Wenger
email: philippe.wenger@ls2n.fr

Kinetostatic analysis and actuation strategy of a planar tensegrity 2-X manipulator

HAL is

Introduction

A tensegrity structure is an assembly of compressive elements (bars) and tensile elements (cables, springs) held together in equilibrium [START_REF] Buckminster | Tensile-integrity structures[END_REF][START_REF] Motro | Tensegrity systems: the state of the art[END_REF]. Tensegrity is known in architecture and art for more than a century [START_REF] Skelton | Tensegrity systems[END_REF] and is suitable for modeling living organisms [START_REF] Levin | The tensegrity-truss as a model for spine mechanics: biotensegrity[END_REF]. Tensegrity mechanisms have been more recently studied for their promising properties in robotics such as low inertia, natural compliance and deployability [START_REF] Arsenault | Kinematic, static and dynamic analysis of a planar 2-dof tensegrity mechanism[END_REF][START_REF] Crane | Kinematic analysis of a planar tensegrity mechanism with pre-stressed springs[END_REF][START_REF] Wenger | Kinetostatic analysis and solution classification of a planar tensegrity mechanism[END_REF]. They are also interesting candidates to design locomotion systems [START_REF] Rieffel | Adaptive and resilient soft tensegrity robots[END_REF][START_REF] Böhm | Spherical mobile robot based on a tensegrity structure with curved compressed members[END_REF][START_REF] Vespignani | Design of superball v2, a compliant tensegrity robot for absorbing large impacts[END_REF][START_REF] Bakker | Design of an environmentally interactive continuum manipulator[END_REF]. A tensegrity mechanism is obtained when one or several elements are actuated. This work falls withing the context of the AVINECK project involving biologists and roboticists with the main objective to model and design bird necks. Birds are capable of using their neck for dextrous tasks interacting with the environment (e.g. a vulture picking food from a carcass), as well as for tasks demanding high force transmissions and accelerations (e.g. the woodpecker hitting a tree trunk). Thus, bird necks are an interesting source of bio-inspiration for designing new manipulators with enhanced performances. Accordingly, a class of planar tensegrity manipulators made of a series assembly of several X-mechanisms i.e. crossed four-bar mechanisms with springs along their lateral sides, has been chosen as a suitable candidate for a preliminary planar model of a bird neck, see Fig. 1. These mechanisms are inspired from the Snelson's X-shape mechanisms [START_REF] Snelson | Continuous tension, discontinuous compression structures[END_REF]. Although simplified because it is planar, this model goes beyond the only available bird neck model in the literature that uses a simple planar articulated linkage [START_REF] Zweers | Motor organization of the avian head-neck system[END_REF], as it can be more easily actuated with cables that play the role of tendons and muscles. Moreover, the center of rotation of the X-mechanism is not fixed, which is the case in most biological joints. Snelson's X-shape mechanisms have been studied by a number of researchers, either as single mechanisms [START_REF] Arsenault | Kinematic, static and dynamic analysis of a planar 2-dof tensegrity mechanism[END_REF][START_REF] Wenger | Kinetostatic analysis and solution classification of a planar tensegrity mechanism[END_REF][START_REF] Boehler | Definition and computation of tensegrity mechanism workspace[END_REF][START_REF] Van Riesen | Dynamic analysis and control of an antagonistically actuated tensegrity mechanism[END_REF] or assembled in series [START_REF] Bakker | Design of an environmentally interactive continuum manipulator[END_REF][START_REF] Aldrich | Time-energy optimal control of hyper-actuated mechanical systems with geometric path constraints[END_REF][START_REF] Arsenault | Kinematic and static analysis of a planar modular 2-dof tensegrity mechanism[END_REF][START_REF] Chen | Analytical computation of the actuator and cartesian workspace boundaries for a planar 2-degree-of-freedom translational tensegrity mechanism[END_REF][START_REF] Furet | Optimal design of tensegrity mechanisms used in a bird neck model[END_REF][START_REF] Furet | Kinematic analysis of planar tensegrity 2-x manipulators[END_REF][START_REF] Furet | Workspace and cuspidality analysis of a 2-x planar manipulator[END_REF]. In this paper, the manipulator may be subject to gravity unlike in [START_REF] Bakker | Design of an environmentally interactive continuum manipulator[END_REF], where the mechanism was used in a snake-like manipulator moving on the ground. This paper is twofold. In the first part, the full kinematic analysis of the manipulator is performed. Since each X-mechanism has a variable instantaneous center of rotation, the symbolic calculation of the inverse geometric problem is more challenging: polynomials of excessive degree may appear if the problem is not Fig. 1: An X-mechanism with two lateral springs (a) and a series assembly mimicking a bird neck (b) tackled with care. A practical method is provided that makes it possible to come up with a quartic polynomial, which turns out to be of minimal degree. The singularities of the manipulator are shown to divide the joint space into only two singularity-free components or aspects. Depending on the values of its geometric parameters, the manipulator at hand may be cuspidal, i.e. it can move from one inverse kinematic solution to another without meeting a singularity. A necessary and sufficient condition is stated for the manipulator to be cuspidal. The effect of joint limits is analyzed and the accessibility regions are further classified according to the admissible configurations of each X-mechanism in these regions. The second part of this paper is dedicated to the static analysis and the determination of the wrench feasible workspace for two different actuation strategies. In the first actuation strategy, the manipulator is actuated with only two symmetric cables. With this actuation strategy, it is shown that the wrench feasible workspace is limited and stiffness cannot be controlled. The second actuation strategy proposed is an over-actuated antagonist strategy that produces a much larger wrench feasible workspace. Moreover, this actuation strategy allows tuning the manipulator stiffness, which can be of interest for safe interaction and collaborative robotics.

Kinematic analysis 2.1 Manipulator description and parametrization

The manipulator consists of a series assembly of two X-mechanisms as shown in Fig. 2. Unlike the original Snelson's X-shape tensegrity that has only two compressive elements (the two crossing bars) and four tensile elements (the four sides), the X-mechanisms used in this paper have four bars in total. The two additional bars replace the tensile elements on the base and top sides. Accordingly, the X-mechanisms used in this study are class-3 tensegrities [START_REF] Skelton | Tensegrity systems[END_REF]. Each X-mechanism i has a base bar and a top bar of length b and two crossed bars of length L, thus defining an antiparallelogram. The mechanism assembly condition is b ≤ L and the particular case b = L defines a degenerate manipulator as shown further. The springs have no impact on the kinematic analysis and have been hidden in Fig. 2 for more simplicity. We need to define a suitable joint variable that describes each mechanism configuration without any ambiguity. In [START_REF] Bakker | Design of an environmentally interactive continuum manipulator[END_REF], the mechanism configuration was chosen as the orientation α i of the upper bar with respect to the base bar, but such a choice is allowed for mechanism motions restricted to -π < α i < π. For a matter of completeness in the kinematic analysis, we are interested in allowing each mechanism to move within its full range beyond its two singularities where it assumes a flat configuration. Since the upper bar makes two turns for one mechanism full turn, α i is not appropriate to define the mechanism configuration. Let us introduce a line segment of length l i that links the middle points of the top and base bars of each mechanism i (shown in dotted line in Fig. 2). The angle between this line and the direction orthogonal to the base bar is referred to as θ i . For the symmetric design considered here, it can be shown that α i = 2θ i . Thus, when the mechanism makes a full turn: -π < θ i < π. Assuming that each mechanism remains always in the crossed-bar assembly mode, the manipulator configuration can thus be unambiguously defined with (θ 1 , θ 2 ). Besides, additional angles, namely φ i and ψ i are considered as they will be useful to derive the inverse kinematics. Since the mechanism sides define an isosceles trapezoid, the length l i of the line segment that links the middle points of the top and base bars can be expressed as follows:

l i (θ i ) = L 2 -b 2 cos 2 (θ i ) (1)
FMANU-JMR-1083, FURET

α 1 1 ¡ 2 ψ 1 ψ 2 l 1 l 2 (x,y) £ 2 ¢ 1 (0,0) ¤ 2 ¤ 1 L b Fig. 2: Manipulator description

Direct Kinematics

To establish the direct kinematic equations, we start with the base frame centered at the middle point of the base bar with the x-axis aligned along this bar. The reference point (x, y) is chosen as the middle point of the top-bar (Fig. 2). The direct kinematic equations of a manipulator with n mechanisms in series can be written as a virtual serial nR manipulator with (varying) link lengths l i and joint angles θ i :

   x = -∑ n i=1 l i (θ i ) sin -θ i + ∑ i j=1 2θ j y = ∑ n i=1 l i (θ i ) cos -θ i + ∑ i j=1 2θ j (2)
For a manipulator made of two mechanisms in series with input variables θ 1 and θ 2 , the direct kinematic equations simplifies to:

x = -l 1 (θ 1 ) sin(θ 1 ) -l 2 (θ 2 ) sin(2θ 1 + θ 2 ) y = l 1 (θ 1 ) cos(θ 1 ) + l 2 (θ 2 ) cos(2θ 1 + θ 2 ) (3)
where l 1 and l 2 were defined in Eq. (1).

Inverse Kinematics

Deriving the so-called characteristic polynomial establishes a univariate polynomial to solve the inverse kinematics. The direct kinematic equations (3) have been derived in their most compact form, namely, the two output variables x and y are expressed as functions of the two input variables θ 1 and θ 2 like in a 2R serial manipulator. However, starting from these equations to derive a characteristic polynomial is not appropriate. Indeed, the square roots appearing in l 1 and l 2 would have to be first cleared out, with the consequence of raising the degree of the resulting characteristic polynomial and providing spurious solutions. Instead, it is better to write x and y as functions of the crossed-bar angles φ 1 and ψ 2 as follows:

x = -b 2 + L cos(φ 1 ) + L cos(ψ 2 + α 1 ) + 1 2 b cos(α 1 + α 2 ) y = L sin(φ 1 ) + L sin(ψ 2 + α 1 ) + 1 2 b sin(α 1 + α 2 ) (4)
In addition, by deriving the anti-parallelogram loop-closure for mechanism i projected on x and y axis, we have:

b + L cos(ψ i ) + b cos(α i ) -L cos(φ i ) = 0 L sin(ψ i ) + b sin(α i ) -L sin(φ i ) = 0 (5)
FMANU-JMR-1083, FURET By rearranging Eqs. [START_REF] Arsenault | Kinematic, static and dynamic analysis of a planar 2-dof tensegrity mechanism[END_REF] for each mechanism, ψ 1 and φ 2 can be eliminated. Finally, the following two loop-closure equations below are obtained:

-2Lb sin(α 1 ) sin(φ 1 ) -2Lb(cos(α 1 ) + 1) cos(φ 1 ) + 2b 2 (cos(α 1 ) + 1) = 0 2Lb sin(α 2 ) sin(ψ 2 ) + 2Lb(cos(α 2 ) + 1) cos(ψ 2 ) + 2b 2 (cos(α 2 ) + 1) = 0 (6)

A set of four equations in four unknowns is thus available. A univariate polynomial can be obtained upon elimination of three of the four unknowns. There are several ways to do so with a computer algebra software. We have used Maple and its Siropa library [START_REF] Jha | Workspace, joint space and singularities of a family of delta-like robot[END_REF]. This library was developed in the frame of a collaborative project on the kinematics of parallel manipulators, see e.g. [START_REF] Manubens | Cusp points in the parameter space of degenerate 3-rpr planar parallel manipulators[END_REF][START_REF] Moroz | On the determination of cusp points of 3-rpr parallel manipulators[END_REF] for more details on the implementation and use of these tools. It contains specific macrofunctions that use efficient algebraic tools such as Groebner bases. Specifically, the Projection function was used to project the system of four equations in order to obtain one single equation in one single variable, chosen here as φ 1 . The half-tangent substitution yields a factored polynomial and one of the factors defines the characteristic polynomial. The characteristic polynomial obtained turns out to be of degree 4 in t=tan(φ 1 /2) and can be written as follows:

a 4 t 4 + a 3 t 3 + a 2 t 2 + a 1 t + a 0 = 0 (7)
where :

a 4 = (b + 1) 2 (b 2 y 2 + x 4 + 2x 2 y 2 + y 4 + 4x 3 + 4xy 2 + 5x 2 + y 2 + 2x) (8a)
a 3 = 4y(b + 1)(2b 2 x + b 2 -2x 2 -2y 2 -4x -1) (8b)
a 2 = 2(b 4 y 2 + b 2 x 4 + 2b 2 x 2 y 2 + b 2 y 4 + b 2 x 2 -10b 2 y 2 + x 4 + 2x 2 y 2 + y 4 -3x 2 + 9y 2 ) (8c) a 1 = 4y(b -1)(2b 2 x -b 2 + 2x 2 + 2y 2 -4x + 1) (8d)
a 0 = (b -1) 2 (b 2 y 2 + x 4 + 2x 2 y 2 + y 4 -4x 3 -4xy 2 + 5x 2 + y 2 -2x) (8e)
Note that L was set equal to 1 for to simplify the calculations without loss of generality. For each solution φ 1 , α 1 is then solved from the first equation in [START_REF] Crane | Kinematic analysis of a planar tensegrity mechanism with pre-stressed springs[END_REF]. Disregarding the solution α 1 = π corresponding to the parallelogram assembly, one obtains:

tan(α 1 /2) = -cos(φ 1 ) + b sin(φ 1 ) (9)
Finally, ψ 2 and α 2 are solved from system (4), which gives two solutions, one of which never satisfies the closed-loop equations [START_REF] Crane | Kinematic analysis of a planar tensegrity mechanism with pre-stressed springs[END_REF]. In total, thus, one comes up with a maximum number of four solutions. The inverse kinematics is solved for a manipulator defined by L = 1 and b = 2/3 at x = 0.03 and y = 1.6. Four solutions are found (see Fig. 3), showing that the characteristic polynomial (7) is of minimal degree.

Singularity analysis and aspects

The Jacobian matrix J of the manipulator, defined such that [ ẋ, ẏ] = J[ θ1 , θ2 ] , is derived from system (3): [START_REF] Vespignani | Design of superball v2, a compliant tensegrity robot for absorbing large impacts[END_REF] with :

J = ∂x ∂θ 1 ∂x ∂θ 2 ∂y ∂θ 1 ∂y ∂θ 2
                   ∂x ∂θ 1 = - -2b 2 cos(θ 1 ) 3 +L 2 cos(θ 1 )+b 2 cos(θ 1 )+2 √ L 2 -b 2 cos(θ 2 ) 2 cos(2θ 1 +θ2) √ L 2 -b 2 cos(θ 1 ) 2 √ L 2 -b 2 cos(θ 1 ) 2 ∂x ∂θ 2 = --sin(2θ 1 +θ 2 )b 2 cos(θ 2 ) sin(θ 2 )-cos(2θ 1 +θ 2 ) cos(θ 2 ) 2 b 2 +L 2 cos(2θ 1 +θ 2 ) √ L 2 -b 2 cos(θ 2 ) 2 ∂y ∂θ 1 = - -2b 2 cos(θ 1 ) 2 sin(θ 1 )+L 2 sin(θ 1 )+2 √ L 2 -b 2 cos(θ 2 ) 2 sin(2θ 1 +θ2) √ L 2 -b 2 cos(θ 1 ) 2 √ L 2 -b 2 cos(θ 1 ) 2 ∂y ∂θ 2 = --cos(2θ 1 +θ 2 )b 2 cos(θ 2 ) sin(θ 2 )-sin(2θ 1 +θ 2 ) cos(θ 2 ) 2 b 2 +L 2 sin(2θ 1 +θ 2 ) √ L 2 -b 2 cos(θ 2 ) 2 (11)
Its determinant is derived and, after clearing the numerator L 2b 2 cos(θ 1 ) 2 L 2b 2 cos(θ 2 ) 2 that never vanishes since b < L, the following expression is obtained:

Det(J) = L 2 -b 2 cos (θ1) sin (θ2) L 2 -2 b 2 cos (θ2) 2 + sin (θ1) cos (θ2) L 2 L 2 + b 2 -2 L 2 b 2 cos (θ2) 2 (12) -2 sin (θ2) cos (θ2) -b 2 cos (θ1) 2 + L 2 L 2 -b 2 cos (θ2) 2 b 2
The solutions to Det(J) = 0 should account for the singularity conditions occuring when the instantaneous centers of rotation (ICR) of the two X-mechanisms are aligned with the end-effector control point. A particular situation arises when the two X-mechanisms are coincident. This situation is equivalent to the fully-folded back configuration of a planar 2R serial manipulator with equal link lengths. When the manipulator is in this singularity, the reference point remains at the origin x = 0 and y = 0. This particular singularity is defined by θ 1 = -θ 2 ± π and it can be verified that Det(J) vanishes in these configurations. The general singularity is like the fully outstretched configuration of a planar 2R serial manipulator. In this singularity, the manipulator cannot produce an instantaneous motion along a line passing through the two ICR and the endeffector control point. Unexpectedly, it was possible to derive an explicit form of this general solution. When cos(θ 2 ) = 0, this solution can be written as: When cos(θ 2 ) = 0, then cos(θ 1 ) = 0, which means that both mechanisms are in their flat singularity. Finally, it can be shown that the second factor of the denominator in (13) never vanishes.

tan(θ 1 ) = -sin(θ 2 ) (2b 2 cos(θ 2 ) 2 -L 2 ) 2 -L 2 b 2 cos(θ 2 ) ((2b 2 cos(θ 2 ) 2 -L 2 ) 2 + 3L 2 b 2 -4 cos(θ 2 ) 2 b 4 ) (13)
The above singularities divide the joint space into singularity-free domains called aspects [START_REF] Borrel | A study of multiple manipulator inverse kinematic solutions with applications to trajectory planning and workspace determination[END_REF]. The singularity θ 1 = -θ 2 ± π produces two lines while the singularity defined by Eq. ( 13) gives a curve, as shown in Fig. 4. In the absence of joint limits, the opposite sides of the square are in fact coincident and the singularity curves divide the joint space into only two aspects. Since the manipulator admits up to four solutions, there are necessarily more than one solution in one of the aspects. In fact, there are two solutions in each aspect, which can be verified by calculating the inverse solutions at any point in a four-solution region. This means that the manipulator can move from one inverse kinematic solution to another without meeting a singularity, namely, it is cuspidal [START_REF] El Omri | How to recognize simply a non-singular posture changing 3-dof manipulator[END_REF][START_REF] Wenger | Cuspidal and noncuspidal robot manipulators[END_REF].

Workspace analysis

When mapped to Cartesian space, the singularity curves define the external and internal boundaries of the workspace. The equation of these curves can be obtained by deriving the discriminant of the characteristic polynomial [START_REF] Wenger | Kinetostatic analysis and solution classification of a planar tensegrity mechanism[END_REF]. By doing so, a polynomial of degree 16 in x and y is obtained. Its expression is quite large and shall not be reported here but can be found in [START_REF]Derivation of a polynomial equation for the boundaries of 2-x manipulators[END_REF]. Figure 5 shows the plot of these singularity curves for three examples. In the middle and right figures, they divide the workspace into three regions. The number of solutions in each regions can be established by counting the number of solutions at any point in each region. In the largest one, the manipulator admits two inverse kinematic solutions. In the two smaller regions, there are four solutions. In the left figure, no internal boundaries exist and the manipulator admits two inverse solutions throughout its workspace. There are two degenerate designs, namely, b = L and b = 0. When b = L, each mechanism degenerates to a bar (the top and base bars coincide and so do their left (resp. right) joints). Thus, these two bars can move with a pure rotation about any of their two revolute joints. Accordingly, the manipulator behaves like two 2R planar linkages of link lengths b and b/2, respectively, with their base joints located at x = -b/2 and x = b/2, respectively (the reason for having b/2 is that the base frame and the end-effector frame are located in the middle of the base and top bars, respectively). Thus, the workspace is the union of two hollow disks of outer (resp. inner) radius 3b/2 (resp. b/2). There are four solutions in the overlapping regions of the non-empty areas of the two hollow disks and two solutions anywhere else (see Fig. 6, left). When b = 0, the manipulator degenerates to a planar 2R linkage and its workspace is a full disk of radius 2L, with two inverse kinematic solutions everywhere, see Fig. 6, right. These two degenerate situations were mentioned for the sake of completeness in the kinematic analysis but they have no practical interest and should be avoided.

The workspace boundaries are associated with the singularities of the 2-X manipulator as a whole. On the other hand, each X-mechanism has its own singularities occurring when it is completely flat. These singularities are defined by θ i = ± π 2 . The role of the X-mechanism singularities is to distinguish mechanism solutions in right side up configurations (the top bar is above the base bar) from those in upside down configurations (the top bar is below the base bar). These two configurations are referred to as configuration up and configuration down, respectively, and illustrated in Fig. 7.

The curves associated with the singularities of the X-mechanisms are now plotted in the workspace, see Fig. 8. The additional curves divide the 2-solution and the 4-solution regions into smaller regions associated with different combinations of mechanism configurations as explained in Fig. 8. For example, the colored region with both | | and // lines means that FMANU-JMR-1083, FURET In the context of a tensegrity manipulator driven with cables along the sides of the mechanisms, the reachable workspace is reduced to those regions where at least one solution exists with both mechanisms up. Those regions have vertical lines in Fig. 8. This point will be discussed in section 3.

Cuspidality condition

As shown above, the workspace may have four-solutions regions or not depending on the values of its geometric parameters. When two 4-solution regions exist, their boundaries have three singular points: a node and two cusps as shown in Fig. 5, right. The existence of cusps shows that the manipulator is cuspidal as a non-singular solution changing motion can be defined by encircling one of the cusps [START_REF] El Omri | How to recognize simply a non-singular posture changing 3-dof manipulator[END_REF][START_REF] Wenger | Cuspidal and noncuspidal robot manipulators[END_REF].

From Fig. 5, the 4-solution regions are smaller when b is smaller. In fact, the four-solution regions may even vanish (Fig. 5, left) and so do also the cusp points. The transition, called swallowtail bifurcation [START_REF] Thomas | On the topological characterization of robot singularity loci. a catastrophe-theoretic approach[END_REF], occurs when the two cusps and the node coincide and define a point with four coincident inverse kinematic solutions. Thus, a geometric condition for having two cusps and four inverse solutions can be established by deriving a condition on the geometric parameters for the 4th-order characteristic polynomial to admit a quadruple root. The latter condition can be stated as follows [START_REF] Wenger | A dh-parameter based condition for 3r orthogonal manipulators to have four distinct inverse kinematic solutions[END_REF]: FMANU-JMR-1083, FURET

     a 4 a 0 -1 4 a 3 a 1 + 1 12 a 2 2 = 0 1 6 a 4 a 2 a 0 -1 16 a 4 a 2 1 -1 16 a 2 3 a 0 + 1 48 a 3 a 2 a 1 -1 216 a 3 2 = 0 1 6 a 4 a 2 -1 16 a 2 3 = 0 ( 14
)
The above system can be simplified by observing that the quadruple root should appear at x = 0 due to symmetry. Solving b from the last two equations in Eq. ( 14) yields b = √ 6/6. When b > √ 6/6 , the manipulator is cuspidal and has four inverse solutions, otherwise it has only two solutions and it is non-cuspidal.

This cuspidality condition was derived based on the assumption that cusps could appear only from a swallowtail bifurcation. While other transitions are possible in general (lips and beaks [START_REF] Thomas | On the topological characterization of robot singularity loci. a catastrophe-theoretic approach[END_REF]), this assumption was set on the basis of intensive observations of the singularity plots for varying values of b. For the sake of rigor, the cuspidality condition was double checked by deriving the general existence condition of cusps. The solution set of cusps in x and y was obtained by eliminating t from the characteristic polynomial and its first and second derivatives [START_REF] Wenger | Cuspidal and noncuspidal robot manipulators[END_REF]. The discriminant of the resulting system was then derived, providing a factored expression in b with 6 factors. Only one factor was found to have real roots, namely 6b 2 -1, which matches the previous condition and confirms that the swallowtail is the only possible bifurcation of the singularity locus.

Keeping in mind that L was set to 1 and taking into account the assembly condition, the cuspidality condition is actually b < L < b √ 6. The following result can thus be stated as follows: a planar 2-degree-of-freedom manipulator made of two identical symmetrical X-mechanisms with crossed link lengths L and base link b is cuspidal and has four inverse kinematic solutions if and only if b < L < b √ 6. To know whether a manipulator is cuspidal or not is very useful for trajectory planning and control [START_REF] Wenger | Cuspidal and noncuspidal robot manipulators[END_REF].

Influence of the joint limits

Unlimited joints were assumed thus far, but for a practical implementation, joint limits will have to be considered. To determine the workspace boundaries in the presence of joint limits, one has to plot the lines associated with these joint limits into the Cartesian space, in addition to the singularity curves. Figure 9 shows how the workspaces with unlimited joints depicted in Fig. 5 change when -π/2 < θ 1 < π/2 and -π/2 < θ 2 < π/2. As compared with Fig. 8, the workspace of the middle figure in Fig. 9 is now reduced to the regions with vertical lines, namely, where both mechanisms are in their right side up configuration. This is because the joint limits restrict the mechanisms from turning upside-down. Note that this restriction is necessary in a remote actuation scheme with active cables since this scheme does not allow passing the flat singular configuration. The number of solutions is determined in each region by solving the inverse kinematics for an arbitrary point. As shown in Fig. 9, the main part of the workspace has only one solution. When b = 2/3 (middle), the 4-solution region is limited to the small upper diamond region (the one with only vertical lines in the zoomed-in view of Fig. 8). Hatching Type of solutions

FMANU-JMR-1083, FURET

| | At least 1 sol. M 1 M 2 // At least 1 sol. M 1 M 2 \\ At least 1 sol. M 1 M 2 == At least 1 sol. M 1 M 2 4 solutions

Actuation strategy

We now focus on the actuation of the tensegrity manipulator at hand. Several actuation schemes are possible, including embedded and remote actuation. Here a remote actuation scheme with active cables is considered. This actuation scheme allows for a light weight design since the actuators are removed from the moving elements. Then, it is necessary to define an approriate number of motors and a suitable way of routing the cables from the motors to the mechanisms. Regarding the number of motors, the minimum number is two since this particular planar tensegrity manipulator has two degrees of freedom. On the other hand, a redundant actuation scheme with four motors would allow each X-mechanism to be antagonistically actuated by active cables along its two sides. With such an antagonistic actuation scheme, we already know that the stiffness of a single mechanism can be tuned and a larger wrench feasible workspace can be obtained [START_REF] Van Riesen | Dynamic analysis and control of an antagonistically actuated tensegrity mechanism[END_REF]. Regarding the routing strategy of the cables, several schemes can be considered, including strut-routed [START_REF] Aldrich | Time-energy optimal control of hyper-actuated mechanical systems with geometric path constraints[END_REF][START_REF] Moored | Analytical predictions, optimization, and design of a tensegrity-based artificial pectoral fin[END_REF] and side-routed [START_REF] Bakker | Design of an environmentally interactive continuum manipulator[END_REF][START_REF] Furet | Optimal design of tensegrity mechanisms used in a bird neck model[END_REF]. Two actuation strategies are studied in this section, one with two actuators and another with four actuators. Since the manipulator is remotely actuated with active cables, each mechanism must be constrained to operate in its up configuration. Accordingly, the rotations range of each mechanism is now restricted to -π 2 < θ i < π 2 .

Static and stiffness modeling

Several hypotheses are considered: (i) the cables are assumed mass-less and infinitely stiff, (ii) all the joints are assumed friction-less and with no clearances, (iii) the cable are assumed to move freely through the nodes without friction and, (iv) all FMANU-JMR-1083, FURET the nonempty part is depicted with vertical lines. The number of solutions is defined by the filling color: none=1, grey=2, red=3 and blue=4 springs are linear tension springs. The static equations relate the input forces and the manipulator configuration. To obtain them, the static equilibrium condition must be established. A convenient way is to derive the total potential energy U total of the tensegrity system at hand. In the absence of external forces and assuming that the cables and the bars are infinitely stiff, U total is composed of the potential energy associated with the springs and the gravity, plus the work produced by the cables, denoted U springs , U gravity and U cables , respectively:

U total = U cables +U springs +U gravity (15)
Let ( f ri , f li ), (l ri , l li ), (k ri , k li ) and (l r0i , l l0i ) refer to the cable forces, cable lengths, spring stiffnesses and spring free lengths of mechanism i, respectively (where subscripts r and l denote left and right, respectively) and let m i and y i denote the mass and the altitude of the gravity center of mechanism i, respectively. Then:

U springs = 1 2 (k r1 (l r1 -l r10 ) 2 + k l1 (l l1 -l l10 ) 2 + k r2 (l r2 -l r20 ) 2 + k l2 (l l2 -l l20 ) 2 ) (16a) U gravity = g(m 1 y 1 + m 2 y 2 ) (16b) U cables = f ( f li , f ri , l li , l ri ) (16c)
The derivation of U cables as a function of the forces and the cable lengths depends on the actuation strategy. It will be detailed in the following sections.

The static equilibrium configurations are the minima of the potential energy. Thus, differentiating Eq. ( 15) with respect to θ 1 and θ 2 yields the relations between the forces in the cables and the configuration of the manipulator:

∂U total ∂θ 1 = 0 ∂U total ∂θ 2 = 0 ( 17
)
The Hessian matrix of the potential energy can be further calculated to analyze the stability of the equilibrium configurations and to define the stiffness matrix K θ of the manipulator [START_REF] Quennouelle | Stiffness matrix of compliant parallel mechanisms[END_REF]:

K θ = K θ 1 K θ 1 θ 2 K θ 1 θ 2 K θ 2 =   ∂ 2 U total ∂θ 2 1 ∂ 2 U total ∂θ 1 ∂θ 2 ∂ 2 U total ∂θ 1 ∂θ 2 ∂ 2 U total ∂θ 2 2   (18)
FMANU-JMR-1083, FURET Fig. 10: Prototype in its rest configuration

The Cartesian stiffness matrix K c of the mechanism is then obtained with the Jacobian matrix J as:

K c = J -K θ J -1 (19)
Controlling the joint or Cartesian stiffness can be of interest in some situations, especially for robots that may interact with obstacles and humans. Consider a task where the robot has to move among obstacles including human workers to a working point where a drilling operation must be carried out. A low joint stiffness during motion is desirable to reduce the impact of collisions, while at the goal point a high Cartesian stiffness in the drilling direction is necessary.

Prototype data

A first prototype is under construction for future experimental studies. It is shown here without its actuation units to provide some practical information, see Fig. 10. The assembly was arranged in space to avoid self collisions. The ratio b/L was mainly decided by considerations regarding the springs extension. Since the minimal and maximal side lengths of each mechanism are Lb and L + b, respectively, b should not be too large to limit the ratio between the minimal and maximal spring lengths. A ratio b/L = 1/2 allows a large choice of springs. Moreover, a compact model that can be easily handled is desired. Accordingly, the bar dimensions were fixed to L = 0.1 m and b = 0.05 m. The workspace along with the number of inverse kinematic solutions in each region are shown in Fig. 11. Note that with these dimensions, the manipulator is cuspidal, but the region where nonsingular posture changing motions are feasible is restricted in the small swallowtail region on top of the workspace. The spring stiffnesses are k l1 = k r1 = 150 N/m and k l2 = k r2 = 90 N/m and their zero-lengths are l l10 = l r10 = 0.047 m and l l20 = l r20 = 0.037 m. The chosen zero-length values insure that the springs are always in tension. The bars were 3D-printed in ABS with a rectangular section of dimension 0.01 m × 0.005 m. Their masses are 0.009 kg for the crossed bars and 0.006 kg for the upper and base bars. All the aforementioned data but the spring stiffnesses shall be used for the calculation of the wrench feasible workspaces in the next sections.

Actuation strategy with two forces

The first actuation strategy studied uses two symmetric cables routed along each side of the mechanisms, as shown in Fig. 12. The work produced by the cables can be written as follows:

U cables = f l (l l1 + l l2 ) + f r (l r1 + l r2 ) (20)
FMANU-JMR-1083, FURET where f l and f r are the cable forces applied on the left and right cable, respectively. Eq. ( 17) can thus be put in the following form:

-f l ∂l l1 ∂θ 1 -f r ∂l r1 ∂θ 1 = k l1 ∂l l1 ∂θ 1 (l l1 -l l10 ) + k r1 ∂l r1 ∂θ 1 (l r1 -l r10 ) + g(m 1 ∂y 1 ∂θ 1 + m 2 ∂y 2 ∂θ 1 ) -f l ∂l l2 ∂θ 2 -f r ∂l r2 ∂θ 2 = k l2 ∂l l2 ∂θ 2 (l l2 -l l20 ) + k r2 ∂l r2 ∂θ 2 (l r2 -l r20 ) + g(m 2 ∂y 2 ∂θ 2 ) ( 21
)
This system is linear in f l and f r . Accordingly, the inverse static model can be solved easily and there is one unique solution ( f l , f r ) for a prescribed (θ 1 , θ 2 ).

The direct static problem, also known as form finding problem in tensegrity structures [START_REF] Skelton | Tensegrity systems[END_REF], requires solving (θ 1 , θ 2 ) for a prescribed ( f l , f r ), which is much more difficult to handle algebraically. Indeed, when deriving the expressions of all derivative terms in system [START_REF] Furet | Workspace and cuspidality analysis of a 2-x planar manipulator[END_REF], complex expressions appear involving square roots of trigonometric expressions both in numerators and denominators. First numerical experiments have revealed the existence of up to seven equilibrium solutions in the ranges -π/2 ≤ θ 1 ≤ π/2, -π/2 ≤ θ 2 ≤ π/2, of which one or two only are stable. The complete algebraic treatment of the direct static problem turns out to be very challenging and is left to future work.

In order to find the set of points (x, y) that the manipulator can reach with this actuation strategy, the joint space (θ 1 , θ 2 ) is discretized and the corresponding points (x, y) and forces ( f l , f r ) are computed. Since the cables can only pull, the forces must be always positive. A minimal positive value is assumed to ensure that the cables are always in tension. Moreover they are bounded by some maximal values depending on the actuators: ( f min , f min ) < ( f l , f r ) < ( f max , f max ). The set of points (x, y) that satisfy these constraints and in which the manipulator can balance a given set of external wrenches is called the wrench feasible workspace (WFW) [START_REF] Ebert-Uphoff | On the connections between cable-driven robots, parallel manipulators and grasping[END_REF]. Here, the only external wrenches considered are the gravity effects. Figure 13 shows the WFW of the manipulator for the prototype data. The WFW is shown in blue. The grey areas were added to show FMANU-JMR-1083, FURET Fig. 13: The WFW in the joint (left) and Cartesian (right) space when the two mechanisms have different springs Fig. 14: The WFW in the joint (left) and Cartesian (right) space when the two mechanisms have identical springs the solutions satisfying only the constraint f l > f min and f r > f min . Clearly, the WFW is much smaller than the geometric workspace. Most points of the WFW lie in the bottom of the workspace, even when the maximum forces are not bounded. The non-symmetry of the WFW comes from the non-symmetry in the spring stiffnesses. Figure 14 shows that when all spring stiffnesses are equal, the size of the WFW is still reduced. This is because in this case the two mechanisms are identical and their behaviors differ only because of the gravity effects. Since these effects are small, most reachable configurations are located around θ 1 = θ 2 . Note that in the absence of gravity and with identical springs, the two mechanisms would have exactly the same behavior and the manipulator would be able to reach only configurations θ 1 = θ 2 , namely, it would have only one degree of freedom.

Actuation strategy with four forces

With the above actuation strategy, the WFW is very small. Moreover, the stiffness cannot be controlled independently of the configuration since there are only two active cables. For a control of both the configuration and the stiffness of the manipulator, four active cables are necessary. Accordingly, a second actuation strategy is studied with two more actuators FMANU-JMR-1083, FURET

f l1 f r1 l l1 ,k l1 l r1 ,k r1 l l2 ,k l2 l r2 ,k r2 (x,y) f r2 f l2
Fig. 15: Second actuation scheme with four cables. The cables driving the second mechanism go along the crossing bars of the first mechanism and then pass along the sides of the second mechanism fixed on the base frame. These two additional motors actuate the first mechanism with two active cables routed along its sides and attached to its upper bar. The second mechanism can be driven with two cables routed along the sides of the two mechanisms like in the first actuation strategy. Another choice is a strut-routed scheme, namely, the cables pass along the crossing bars of the first mechanism. Both solutions have their own benefits and drawbacks [START_REF] Moored | Analytical predictions, optimization, and design of a tensegrity-based artificial pectoral fin[END_REF]. One important drawback of the side-routed scheme is that it creates couplings in the cable lengths. Moreover, it tends to increase the stiffness of the first mechanism, since it is affected by the four forces. Accordingly, a strut-routed scheme is chosen (see Fig. 15). For this actuation strategy, The work produced by the cables is:

U cables = f l1 l l1 + f r1 l r1 + f l2 l l2 + f r2 l r2 (22)
As there are four forces for only two degrees of freedom, system of Eqs. ( 17) has twice-infinitely many solutions. The force redundancy can be solved by imposing a minimal or a maximal stiffness value of two elements in K θ , depending on whether safe interactions are desirable or a high stiffness is required. Since the forces appear only in the diagonal elements of K θ (the non-diagonal terms depend only on the gravity), the two diagonal terms K θ 1 and K θ 2 can be chosen. To compute the WFW, the joint space is scanned for each configuration (θ 1 , θ 2 ), and the set of forces ( f l1 , f l2 , f r1 , f r2 ) with minimal and maximal values of K θ 1 and K θ 2 is computed. Since each mechanism i is actuated with two antagonistic forces ( f li , f ri ) independently from the other mechanism, the minimal stiffness K θ i is obtained when one of the two forces ( f li , f ri ) vanishes. An algorithm was developed to determine the two zero forces associated to (θ 1 , θ 2 ) and to compute the other two forces using Eq. [START_REF] Arsenault | Kinematic and static analysis of a planar modular 2-dof tensegrity mechanism[END_REF]. If the computed forces fall within their bounds, the corresponding stiffnesses K θ 1 and K θ 2 are then calculated using Eq. [START_REF] Chen | Analytical computation of the actuator and cartesian workspace boundaries for a planar 2-degree-of-freedom translational tensegrity mechanism[END_REF]. A similar algorithm is used to compute the maximal stiffnesses, but this time at least two forces are set to their maximal value. Figure 16 (resp. Fig. 17) shows the wrench-feasible configurations in the joint space along with the minimal (resp. maximal) joint stiffnesses. These figures show that the joint stiffness of the i th mechanism mainly depends on its own θ i position, but there is also a small influence of the other mechanism, which is due to the gravity effects. It can be seen here that all the configurations (θ 1 , θ 2 ) are inside the WFW. Some configurations near the joint limits have negative stiffness, which means that the manipulator will be in an unstable equilibrium there. This is due to the closeness to the flat singularities, that are unstable. Such configurations could be avoided by reducing the joint limits.

For each (θ 1 , θ 2 ), the corresponding point (x, y) is calculated with Eq. ( 3) and the Cartesian stiffness K c is computed using Eq. [START_REF] Furet | Optimal design of tensegrity mechanisms used in a bird neck model[END_REF]. The maximal and minimal values of the diagonal elements K x and K y , which define the stiffness in the x and y direction, respectively, are computed like above (it can be shown from their algebraic expressions that K x and K y attain their maximum and minimum values for the same forces as K θ 1 and K θ 2 ). Figure 18 (resp. Fig. 19) shows the WFW along with the minimal and maximal Cartesian stiffnesses in the x and y directions. The Cartesian stiffness is naturally influenced by the manipulator singularities as well as by the mechanism singularities. They even become infinite in the direction defined by the two instantaneous centers of rotation and the end-effector point. To avoid infinite values, the calculation of the Cartesian stiffness was discarded around the singularities. It is worth noting that there is generally a ratio of 10 between the minimal and the maximal stiffnesses both in the joint space and in the Cartesian space, which allows tuning the stiffness in a rather large range. Note that the stiffness mapping was conducted in the x and y directions, but any other directions could have been considered.

FMANU-JMR-1083, FURET This paper has reported on a preliminary study of a project involving biologists that aims at producing a bird neck model as a bio-inspired manipulator design with enhanced performances. Accordingly, a 2-degree-of-freedom tensegrity manipulator made of two X-mechanisms in series has been studied. Although simple at a first glance, this manipulator turned out to have some challenging kinematics issues, such as an inverse kinematics that requires an adequate solving methodology and a cuspidal behavior. On the other hand, it lends itself well to remote actuation schemes with cables and springs. The first part was dedicated to a pure kinematic analysis, where a methodology to solve the inverse kinematics was proposed, the singularities were determined and the workspace with and without joint limits was analyzed in details. The manipulator was shown to have two or four solutions depending on its geometric parameter values. When it has four solutions, the manipulator was shown to be cuspidal, a result that was never encountered before in a planar 2-dof manipulator, to the best of the authors' knowledge. A cuspidality condition based on the geometric parameters was established. In the second part, two actuation schemes with remote motors, lateral springs and cables have been proposed. The first scheme uses only two motors, which is enough from a pure kinematic point of view. However, this scheme turned out to produce a limited wrench feasible workspace (WFW). The second actuation strategy uses four remote motors, thus introducing actuation redundancy. In addition to producing a much larger WFW, this redundant actuation strategy allows tuning the stiffness in the workspace. This feature is of high interest for tasks requiring both safe interactions and high forces transmission. Only static analyses have been reported in this paper. The dynamic behavior is under study and first simulation results with a dynamic control FMANU-JMR-1083, FURET law reveal satisfactory responses [START_REF] Fasquelle | Dynamic modeling and control of a tensegrity manipulator mimicking a bird neck[END_REF]. Experiments are now being conducted on a prototype and shall be presented in a future paper. First results show that friction, which was neglected in this study, has signification effects. It can be taken into account by adding a static force as an offset in the active forces. In future work, the planar model will be completed with more mechanisms in series. Kinematic redundancy will allow avoiding obstacles, and even shaping around them, where contacts with obstacles will be modelled with external forces. Finally, a 3D model will be proposed, which will account for the additional degree of freedom existing between two subsequent vertebrae in the real bird neck.

FMANUFig. 3 :

 3 Fig. 3: The four inverse solutions at x = 0.03 and y = 1.6 (L = 1 and b = 2/3)

FMANUFig. 4 :

 4 Fig. 4: Singularity curves in the joint space and aspects when L = 1 and b = 2/3

Fig. 5 :Fig. 6 :

 56 Fig. 5: Workspace boundaries when L = 1 and b = 2/5 (a), b = 2/3 (b), b = 9/10 (c). Four-solutions regions are filled in gray

Fig. 7 :

 7 Fig. 7: X-mechanism in configuration up (a) and down (b)

Fig. 8 :

 8 Fig. 8: workspace regions with different numbers and types of solutions (L = 1, b = 2/3. M i (resp. M i ) means that mechanism i is in configuration down (resp. up)

Fig. 9 :

 9 Fig. 9: workspace boundaries when -π/2 < θ 1 , θ 2 < π/2, and for L = 1, b = 2/5 (left), b = 2/3 (center), b = 9/10 (right).the nonempty part is depicted with vertical lines. The number of solutions is defined by the filling color: none=1, grey=2, red=3 and blue=4

Fig. 11 :Fig. 12 :

 1112 Fig.11: Workspace of the manipulator with the prototype dimensions. The nonempty areas are filled with vertical lines. The number of solutions is defined by the filling color: none=1, grey=2 and blue=4

Fig. 16 :

 16 Fig. 16: Minimal joint stiffnesses

Fig. 18 :Fig. 19 :

 1819 Fig. 18: Wrench feasible workspace with the minimal Cartesian stiffness mapping

2 3 4 5 6

 23456 mechanism with two lateral springs (a) and a series assembly mimicking a bird neck (b) . . . . . . . 2 Manipulator description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 The four inverse solutions at x = 0.03 and y = 1.6 (L = 1 and b = 2/3) . . . . . . . . . . . . . . . . . . . . 5 Singularity curves in the joint space and aspects when L = 1 and b = 2/3 . . . . . . . . . . . . . . . . . . 6 Workspace boundaries when L = 1 and b = 2/5 (a), b = 2/3 (b), b = 9/10 (c). Four-solutions regions are filled in gray . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 Workspace boundaries for the two degenerate designs b = L = 1 (a) and b = 0, L = 1 (b). Four-solutions regions are filled in gray . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 7 X-mechanism in configuration up (a) and down (b) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 8 workspace regions with different numbers and types of solutions (L = 1, b = 2/3. M i (resp. M i ) means that mechanism i is in configuration down (resp. up) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 9 workspace boundaries when -π/2 < θ 1 , θ 2 < π/2, and for L = 1, b = 2/5 (left), b = 2/3 (center), b = 9/10 (right). the nonempty part is depicted with vertical lines. The number of solutions is defined by the filling color: none=1, grey=2, red=3 and blue=4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 10 Prototype in its rest configuration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 11 Workspace of the manipulator with the prototype dimensions. The nonempty areas are filled with vertical lines. The number of solutions is defined by the filling color: none=1, grey=2 and blue=4 . . . . . . . . . . 12 12 First actuation strategy with two cables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 13 The WFW in the joint (left) and Cartesian (right) space when the two mechanisms have different springs . . 13 14 The WFW in the joint (left) and Cartesian (right) space when the two mechanisms have identical springs . . 13 15 Second actuation scheme with four cables. The cables driving the second mechanism go along the crossing bars of the first mechanism and then pass along the sides of the second mechanism . . . . . . . . . . . . . 14 16 Minimal joint stiffnesses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 17 Maximal joint stiffnesses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 18 Wrench feasible workspace with the minimal Cartesian stiffness mapping . . . . . . . . . . . . . . . . . . 16 19 Wrench feasible workspace with the maximal Cartesian stiffness mapping . . . . . . . . . . . . . . . . . . 16 FMANU-JMR-1083, FURET

Acknowledgements

This work was conducted with the support of the French National Research Agency (AVINECK Project ANR-16-CE33-0025).

FMANU-JMR-1083, FURET

Nomenclature

L Diagonal bar length b Bottom and top bar length θ i Median bar orientation of mechanism i l i Median bar length of mechanism i J Jacobian matrix of the manipulator k ri Right spring stiffness of mechanism i k li Left spring stiffness of mechanism i l ri Right spring length of mechanism i l li Left spring length of mechanism i l ri0 Right spring zero free length of mechanism i l li0 Left spring zero free length of mechanism i f ri Force along right spring of mechanism i f li Force along left spring of mechanism i g Gravity acceleration m i Mass of mechanism i y i Height of CoM of mechanism i U total Total potential energy of manipulator K θ Hessian matrix of manipulator K c Cartesian stiffness matrix of manipulator


Sandrine Vaton

Fractal

come

Fractal' versus 'Markov' models of traffic (conférence invitée)


S Bustamante

J Taniguchi

M V V Silva

R C C Petry

Montone

C V Z Cipro

P Bustamante

S Taniguchi

M V Petry

R C Montone

Seabird colonies as relevant sources of

Keywords: Antarctica, secondary sources, stable isotopes, POPs, PCBs, Organochlorine pesticides

-Introduction

As presented in part 1 of our work (Cipro et al., 2018a), "global distillation", the mechanism by which contaminants volatilise in warmer locations with subsequent atmospheric long-range transport and condensation in colder regions is majorly held responsible as the main, or even the only pollutant input process in Polar environments, whilst the biologically mediated transport is frequently not taken into account [START_REF] Blais | Arctic Seabirds Transport Marine-Derived Contaminants[END_REF][START_REF] Blais | Biologically Mediated Transport of Contaminants to Aquatic Systems Critical Review Biologically Mediated Transport of Contaminants to[END_REF]. However, literature presents evidences suggesting the role of seabird colonies as relevant secondary sources at the local or regional level for both organic and inorganic contaminants (Cipro et al., 2018a[START_REF] Cipro | Organic pollutants and their correlation with stable isotopes in vegetation from King George Island, Antarctica[END_REF][START_REF] Roosens | Penguin colonies as secondary sources of contamination with persistent organic pollutants[END_REF]. The particular fact that several species of seabirds assemble in colonies containing a large number of individuals (penguins, for instance, may form colonies of hundreds of thousands of individuals, e.g., [START_REF] Ballance | An energetic correlate between colony size and foraging effort in seabirds, an example of the adélie penguin pygoscelis adeliae[END_REF] excreting on land, notably during the reproductive season, results in a substantial input regarding the terrestrial ecosystems in and around such colonies. This input from seabird colonies added to the organic matter deriving from eggs [START_REF] Brasso | Penguin eggshell membranes reflect homogeneity of mercury in the marine food web surrounding the Antarctic Peninsula[END_REF], bird tissues and prey present in these structures [START_REF] Emslie | Ornithogenic soils and the paleoecology of pygoscelid penguins in antarctica[END_REF] represent the major source of nutrients for Antarctic terrestrial ecosystems [START_REF] Smykla | Zonation of Vegetation Related to Penguin Rookeries on King George Island , Maritime Antarctic[END_REF]. This is verified for both extant and abandoned colonies [START_REF] Liu | δ13C and δ15N in the ornithogenic sediments from the Antarctic maritime as palaeoecological proxies during the past 2000 yr[END_REF][START_REF] Outridge | Seabird transfer of nutrients and trace elements from the north water polynya to land during the midholocene warm period, carey islands, northwest Greenland[END_REF][START_REF] Tatur | Formation of abandoned penguin rookery ecosystems in the maritime Antarctic[END_REF]. The availability of the nutrients resulting from this input induces a pattern of green algae/cyanobacteria, Antarctic hair-grass, mosses and lichens, in this order, as the distance to the colony increases [START_REF] Smykla | Zonation of Vegetation Related to Penguin Rookeries on King George Island , Maritime Antarctic[END_REF]. In this regard, water availability also plays an important role considering the differences in water dependency from these organisms but also seasonal melting and runoff water are known to have a key influence in marine primary production [START_REF] Anderson | Nutrient fluxes from water to land: seabirds affect plant nutrient status on Gulf of California islands[END_REF] and in carrying organic contaminants (Cipro et al., 2017b) to be absorbed by marine life (Cipro et al., 2018b). This introduction of nutrients brings contaminants along, not only by a direct influence of the excretion [START_REF] Rudolph | Persistent organic pollutants and porphyrins biomarkers in penguin faeces from Kopaitic Island and Antarctic Peninsula[END_REF], but also through cycles of resuspension and redeposition, for instance.

The suitability of the matrixes for this objective (the assessment of seabird colonies as secondary contaminants sources) is largely confirmed by literature (e.g. [START_REF] Borghini | Organochlorine pollutants in soils and mosses from Victoria Land (Antarctica)[END_REF][START_REF] Cocks | On the uptake of ornithogenic products by plants on the inland mountains of Dronning Maud Land, Antarctica, using stable isotopes[END_REF][START_REF] Negoita | Distribution of polychlorinated biphenyls (PCBs) and organochlorine pesticides in soils from the East Antarctic coast[END_REF][START_REF] Tarcau | Organochlorine pesticides in soil, moss and tree-bark from North-Eastern Romania[END_REF]. However, differently from the trace elements presented in the first part of the present study (Cipro et al. 2018a), the Persistent Organic Pollutants (POPs) analysed in this second part do not present natural environmental sources, but only anthropic ones. Due to their physico-chemical properties (mainly stability, lipophilicity and vapour pressure), POPs 1) are transported through long distances through the atmospheric pathway; 2) undergo an increase of their concentrations during the lives of the organisms (i.e. bioaccumulation process) and also according to their trophic position (i.e. biomagnification process). The cold conditions of Antarctic environments favour the persistency of these compounds with respect to temperate and tropic environments, particularly because of the seasonal storage and consumption of lipids in these environments [START_REF] Cipro | Persistent organic pollutants in bird, fish and invertebrate samples from King George Island, Antarctica[END_REF][START_REF] Loganathan | Time perspectives of organochlorine contamination in the global environment[END_REF].

Taking all the previous considerations into account, in an analogous manner to the first part of the study, the work had as an aim, to assess seabird colonies as secondary sources of selected POPs using both soil and vegetation as indicators in several Antarctic terrestrial ecosystems considering geographical variation, distance to the colonies and the use of different seabird species with different feeding ecology strategies.

-Material and methods

-Sampling

Since this is a two parts study, the present work was performed with the same samples used in the first one (Cipro et al., 2018a). Samples of lichen (n=55), mosses (n=58) and soils (n=37), were collected in 13 locations within the South Shetlands Archipelago in the austral summers of 2013-14 and 2014-15 during cruises of the Brazilian Antarctic Programme and divided in two types: "colony" (within the colony itself for soil and within the colony or bordering such colony for vegetation) and "control" (at least 150m away from the respectively closest colony. Lichens and soil were collected in higher elevations than the colony, whereas mosses, due to their water-dependency, were collected in lower sites. All samples for POP analyses were collected with steel gear (spoons and tweezers) previously rinsed with n-hexane.

All samples for stable isotope analysis (SIA) were collected with plastic gear (spoons and tweezers) firstly rinsed in an acid bath (35 mL L -1 nitric acid, 50 mL L -1 hydrochloric acid), then further rinsed with Milli-Q water. Soil samples were dug no deeper than 5cm. They were chosen at the field mainly by their availability, avoiding fresh faeces, which could mask the results for the formed ornithogenic soil. Samples were collected in triplicates whenever possible.

In some cases, samples were collected in more than one spot in each area of interest, which was described by the "main location" and "specific location" labels. Thus, the comparison between control and colony sites was made within the same specific location, when available, otherwise the closest control sample within the same main location was adopted.

Once collected, samples for POPs analyses were stored in previously combusted aluminium foil and samples for SIA were stored in hermetically sealed plastic bags. They were frozen aboard (-20°C) and kept frozen until arrival at the Marine Organic Chemistry Laboratory (LabQOM, University of São Paulo/Brazil), where all POPs analyses were performed and, also, where the SIA samples underwent lyophilisation. Vegetation species were identified by Prof.

Jair Putzke (Universidade Federal do Pampa, São Gabriel, Brazil). SIA samples were then sent to the University of La Rochelle, France, where they were ground to a fine powder in a ceramic mortar before analyses. Soil samples were sifted (1mm mesh) after freeze-drying in order to remove larger rock fragments and debris. Soil samples were not ground in any part of the process whatsoever in order to assure that only the granulometry of choice underwent analyses.

-Analyses

Analyses for POPs were performed for three classes of contaminants: polychlorinated biphenyls (PCBs), organochlorine pesticides (OCPs) and polybrominated diphenyl ethers (PBDEs).

Sample extraction followed a previously described method [START_REF] Cipro | Persistent organic pollutants in bird, fish and invertebrate samples from King George Island, Antarctica[END_REF], which resulted from an adaptation from literature [START_REF] Loganathan | Time perspectives of organochlorine contamination in the global environment[END_REF]. Briefly, wet samples (10g for vegetation and colony soil, 20g for control soil) were manually ground with anhydrous Na2SO4 and surrogate (PCB 103) was added before extraction in a Soxhlet apparatus for 8 hours with 80ml of n-hexane and methylene chloride (1:1, v/v). During the extraction, metallic copper was added to the flasks in order to retain sulphur. The extract was then concentrated to 1ml and cleaned up in a column filled (from top to bottom) with 16g alumina and 8g silica gel (both 5% deactivated with water). A second rotoevaporation (up to 900 l) followed and finally, internal standard (100ng of TCMX, used to estimate surrogate recovery) was added to the purified extract prior to injection in the gas chromatograph. The same extract is used for all POPs analyses.

PCBs and PBDEs analyses were performed in an Agilent 6890 Plus attached to the MS 5973N Mass Selective Detector (GC/MS) in SIM (selective ion monitoring) mode, with a HP-5MS column (30mx250 mx0.25 m, internally coated with 5% phenyl-95% dimethylpolysiloxane). Helium was used as carrier gas at a constant flow (1.1 mL min -1 ).

Injection volume was 1 L in splitless mode. Injector, interface and ion source operated at 280°C, 280°C and 300°C respectively. The oven ramp was programmed as follows: 75°C for 3 minutes, then rising at 15°C min -1 up to 150°C, then rising at 2°C min -1 up to 260°C, then rising at 20°C min -1 up to 300°C and remaining as so for 10 minutes for a total run time of 75 minutes.

The analysed PCBs congeners were the IUPAC # 8,18,28,31,33,44,49,52,56/60,66,70,74,77,81,87,95,97,101,105,114,118,123,126,128,132,138,141,149,153,156,157,158,167,169,170,174,177,180,183,187,189,194,195,201,203,206 and 209. In turn, the analysed PBDEs congeners were the IUPAC# 28, 47, 99, 100, 153, 154 and183. OCPs analyses were run in a gas chromatograph equipped with an electron capture detector (GC-ECD, Agilent Technologies, model 6890N). Hydrogen was used as carrier gas at a constant pressure (13.0 psi). The injector operated in splitless mode and was kept at 300°C.

The capillary column used was the same HP-5MS. The detector operated at 320°C using N 2 as makeup gas at a flow rate of 60mL min -1 . The oven was programmed as follows: starting at 60°C, then rising at 5°C min -1 up to 150°C, remaining as so for 6 minutes, then rising at 1°C min -1 up to 200°C, then rising at 18°C min -1 up to 300°C remaining as so until a final run time of 90 minutes. The compounds analysed were the HCB (hexachlorobenzene), DDTs (Dichlorodiphenyltrichloroethane DDT, Dichlorodiphenyldichloroethylene DDE and Dichlorodiphenyldichloroethane DDD in op' and pp' configurations), HCHs (Hexachlorocyclohexane in α, β, γ, and δ isomers), chlordanes (αand γ-chlordane, oxichlordane, heptachlor, and heptachlor epoxide), mirex, metoxychlor, endosulfan and drins (aldrin, dieldrin, isodrin and endrin). OCPs quantifications were performed in GC-ECD only after confirmation in GC/MS, which decreased the number of detected compounds, as further discussed.

Method detection limits (MDLs) were calculated as being three times the standard deviation after seven blank replicates. Spiked matrices were recovered within the acceptance ranges (i.e. 40-130% for at least 80% of the spiked analytes) suggested by QA/QC standards [START_REF] Wade | Use of standards and reference materials in the measurement of chlorinated hydrocarbon residues[END_REF]. Method validation was performed using NIST SRM 1945. Blanks were included in every analytical batch (usually 10-12 samples) and all data were blanksubtracted. Average surrogate recovery was 93.27%.

SIA was performed as described in the first part of this work (Cipro et al., 2018a), with the same restrictions, especially the decarbonation of soil samples. Due to the negligible lipid content and influence in the results verified during previous studies (Cipro et al., 2017a[START_REF] Cipro | Organic pollutants and their correlation with stable isotopes in vegetation from King George Island, Antarctica[END_REF], samples have not undergone delipidation. Soil samples, however, have undergone a decarbonation procedure to avoid the interference of carbonates, which are depleted in 13 C and could represent a bias in data interpretation. The decarbonation was as follows: up to 100mg of ground sample were put in a glass vial, 1mL of HCl 0.1N were added and the vial was placed in a microwave bath, observing bubble formation as an indicator of carbonates digestion. After one minute, 100µl of the same acid was added to verify that no more bubble formation was shown. If that was the case, this first step was repeated many times as needed. Next, the vials were put in a 60°C dry bath (Techne) coupled to an evaporation system consisting of tubes gently blowing filtered analytical quality compressed air into the vials overnight in order to evaporate the liquid acidic phase. The following day, samples were recovered with 1mL of Milli-Q water, homogenised in a microwave bath for 1 minute, frozen and finally lyophilised.

SIA themselves were performed as follows: an aliquot of ground prepared sample (0.8 up to 1.5 mg for vegetation, 1.5 to 8 mg for colony soil, 8 to 16 mg for control soil) was encapsulated in tin cups and injected in a Thermo Scientific Delta V AdvantageConFlo IV interface (NoBlank and SmartEA) coupled to a Thermo ScientificFlash EA1112 Elemental Analyser. Pee Dee Belemnite and atmospheric nitrogen were used as standards for calculation of δ 13 C and δ 15 N, respectively. Based on replicate measurements of internal laboratory standards, experimental precision is of ±0.15‰ and ±0.20‰ for δ 13 C and δ 15 N, respectively.

-Statistics

Statistics were performed using Microsoft Excel, Minitab 18 and Statsoft Statistica 10 and 13. Before analyses, data were checked for normality of distribution and homogeneity of variances using Shapiro-Wilk and Brown-Forsythe tests, respectively. Statistically significant results are set at =0.05 unless stated otherwise. Differences of colony and control datasets were assessed by means of t-tests and Pearson/Spearman correlations are chosen due to their linear/monotonic nature and not necessarily to parametric/non-parametric analyses.

Normal distribution and log-link function generalized linear models (GLM) were built as follows: POPs concentrations as the dependent variable; matrix, main location, specific location, colony species and matrix species as categorical factors; δ 13 C, δ 15 N, %C and %N as continuous predictors. Biologically relevant models were constructed incorporating the different variables and their interactions taking the care of removing continuous variables that were significantly correlated before model building in each dataset and its subsets as well.

Model selection was based on Akaïke's Information Criteria adjusted for small sample sizes (AICc). The model with the lowest AICc value was considered to be the most accurate. Models with AICc values differing by less than 2 have a similar level of support in the data, and the model including the least number of parameters can be regarded as the most accurate, according to the principle of parsimony [START_REF] Burnham | Model Selection and Multimodel Inference[END_REF]. Overall model support was assessed using Akaike weights (wi), following [START_REF] Johnson | Model selection in ecology and evolution[END_REF]. Residual (R 2 adj)

analyses should be restricted to description and not be used in model selection [START_REF] Burnham | Model Selection and Multimodel Inference[END_REF]. Only models with AICc values differing by less than 10 are presented.

-Results and discussion

The most present compounds in both vegetation and soil samples were HCB and PCBs. Nevertheless, the full dataset is present as supplementary material (Table S1).

Regarding the concentrations presented in Table 1, some patterns appear. In a general way, for control samples, lichens and mosses presented similar PCBs concentrations and both were one order of magnitude above those found in control soil. Considering colony samples, the pattern is somewhat similar: lichens presenting concentrations slightly higher than mosses and both one order of magnitude higher than soil in average again. Therefore, the main difference in PCBs distribution between colony and control samples would be a relative increase in lichens concentration in colony samples, meaning that, proportionally, most of the input derived from the colony ends up absorbed by this matrix, possibly by volatilisation and/or resuspension. HCB, in turn, presented mean concentrations in control samples in mosses around 50% higher than those found in lichens; which in turn presented values that were double those found in soil, all of them within the same order of magnitude. In colony samples, however, HCB presented overall concentrations in soil samples one order of magnitude above those found in mosses (with the exception of P. papua from Hannah Point likely due to some local effect such as water percolation); and mosses, on the other hand, presented concentrations one order of magnitude higher than lichens. This leads to the conclusion that firstly, HCB is proportionally more excreted than PCBs. Moreover, HCB dominated seabird profiles are reported by literature [START_REF] Carravieri | Wandering Albatrosses Document Latitudinal Variations in the Transfer of Persistent Organic Pollutants and Mercury to Southern Ocean Predators[END_REF] and, in this regard, similar data from chicks [START_REF] Carravieri | From Antarctica to the subtropics: Contrasted geographical concentrations of selenium, mercury, and persistent organic pollutants in skua chicks (Catharacta spp.)[END_REF] evidentiate this phenomenon since most of their POPs burden comes from what is not accumulated by the mother. In addition to that, several previous works demonstrate proportionally lighter profiles in excreta [START_REF] Roosens | Penguin colonies as secondary sources of contamination with persistent organic pollutants[END_REF][START_REF] Rudolph | Persistent organic pollutants and porphyrins biomarkers in penguin faeces from Kopaitic Island and Antarctic Peninsula[END_REF]. Another consideration is yet to be made: the very high levels found in the P. atriceps soil colony samples is likely due to the fact that they form dense colonies in high and dry spots, so very little water percolation occurs when compared to penguin colonies for instance.

Secondly, volatilisation and/or resuspension play a proportionally lesser role in this case when compared to PCBs, likely due to a higher water solubility (therefore a lower Kow of HCB than most PCBs) of the former regarding the latter and also due to a higher water dependency of mosses compared to lichens (see [START_REF] Cipro | Organic pollutants and their correlation with stable isotopes in vegetation from King George Island, Antarctica[END_REF]. This interpretation, however, hangs on the lack of data concerning the absorption mechanism of POPs in mosses. Simple diffusion can explain our results, but there is probably some non-linearity in the absorption process. The results for the even more water soluble HCHs could enlighten this issue, but unfortunately they have remained under the method detection level.

In order to go deeper in the understanding of POPs occurrence and dynamics, correlation analyses and PCA have been performed. Spearman correlation results are presented in Table 2, whereas PCA is presented as supplementary material (Fig. S1).

Lichen in control sites presented two significant correlations for POPs. Firstly, a positive one of PCBs and %N, meaning that a common source is possibly shared, likely atmospheric deposition from long range transport [START_REF] Montone | PCBs and chlorinated pesticides (DDTs, HCHs and HCB) in the atmosphere of the southwest Atlantic and Antarctic oceans[END_REF][START_REF] Montone | PCBs in the atmosphere of King George Island, Antarctica[END_REF], even if the role of colonies as sources via resuspension/volatilisation cannot be discarded at this point. A negative correlation for HCB and %C was also found for this subgroup, leading to two possible, nonmutually excluding hypotheses: the range of the resuspended/volatilised HCB from sources is greater than that of organic carbon and also that some organic matter decomposition process might play a role as hypothesized in the previous work (Cipro et al., 2018a). This would mean that HCB might be liberated in the atmosphere but the organic carbon might be trapped by some previous step of its biogeochemical cycle. Lichen in colony sites, however, presented several significant correlations. PCBs correlated positively with %C and negatively with δ 15 N whereas HCB presented the same ones and also a negative correlation with %N. The first two ones can be explained by the decomposition of the animal derived organic matter. The kinetic fractionations during the decomposition of deposited urea and uric acid lead to the formation of starkly 15 N-depleted highly volatile ammonia whereas the remaining ammonium, much less volatile, is conversely 15 N-enriched [START_REF] Heaton | Isotopic studies of nitrogen pollution in the hydrosphere and atmosphere: a review[END_REF]. Therefore it results that vegetation around the excrement zone are more exposed to 15 N-enriched inorganic nitrogen whereas upland sites are more exposed to 15 N-depleted nitrogen (Cipro et al., 2018a[START_REF] Cipro | Organic pollutants and their correlation with stable isotopes in vegetation from King George Island, Antarctica[END_REF], as indicated by the results found. Moreover, and as a consequence of the previous statement, animal-derived nitrogen uptake is associated with large δ 15 N ranges, with typical values ranging between 6 and 26‰ (e.g. [START_REF] Cipro | Organic pollutants and their correlation with stable isotopes in vegetation from King George Island, Antarctica[END_REF][START_REF] Cocks | On the uptake of ornithogenic products by plants on the inland mountains of Dronning Maud Land, Antarctica, using stable isotopes[END_REF][START_REF] Erskine | Subantarctic Macquarie Island -a model ecosystem for studying animal-derived nitrogen sources using 15 N natural abundance[END_REF][START_REF] Mizutani | Nitrogen and carbon isotope ratios in seabird rookeries and their ecological implications[END_REF]. The negative correlation between HCB and %N is likely explained in an analogous manner to the previous negative correlation between this contaminant and %C in control sites, i.e. the contaminant might be released to the atmosphere but this time, we hypothesise that the nitrogen might be trapped by some previous step of its biogeochemical cycle. Some of these same patterns were found for SIA and trace elements, particularly the more volatile Hg in the previous study (Cipro et al., 2018a).

Moss in control sites presented no significant correlation for any of the contaminants.

In colony sites, however, a negative correlation between PCBs and %N and a positive one between HCB and δ 13 C were detected. The first one, besides the previously discussed occurrence in lichens, was also found in several of the trace elements (Co, Fe, Mn, Ni, Pb, and V) analysed in the first part of the study (Cipro et al., 2018a). Some of these, particularly Mn and Pb, have already been reported with comparatively higher concentrations in Gentoo penguins (Pygoscelis papua) faeces [START_REF] Celis | Trace Element Concentrations in Biotic Matrices of Gentoo Penguins (Pygoscelis Papua) and Coastal Soils from Different Locations of the Antarctic Peninsula[END_REF] and, therefore, might confirm the role of colonies as secondary sources for PCBs in an analogous manner. The correlation between HCB and δ 13 C is at a first glance, somewhat contrary to the previous literature that suggests lower δ 13 C values are related to moister habitats and marine influence [START_REF] Cipro | Organic pollutants and their correlation with stable isotopes in vegetation from King George Island, Antarctica[END_REF][START_REF] Huiskes | Stable isotope ratios as a tool for assessing changes in carbon and nutrient sources in Antarctic terrestrial ecosystems[END_REF][START_REF] Lee | Carbon and nitrogen isotope composition of vegetation on King George Island, maritime Antarctic[END_REF], since mosses are more water dependant than lichens, negative correlations were to be expected. However, this intra-group correlation means that the samples in dryer locations are more exposed to HCB, which is consistent with its comparatively higher volatility.

Soil presented no significant correlations for any of the contaminants in the present study neither for control nor for colony sites.

PCA results (Fig S1, supplementary material) reinforce the previously discussed correlations on the subsets of samples. Nevertheless, when the broader data is taken into account two interesting trends show up: the pattern for all lichen samples taken together resembles closely that for its control sites, whereas the pattern for all moss samples taken together resembles even more closely that for its colony sites. This likely means that the mechanisms for POPs exposure in control sites are more representative for the entirety of the lichens, i.e., long-range atmospheric transport is probably the major force in this case. Mosses, on the other hand, that presented patterns for the whole dataset even closer to the ones regarding colony samples, allow an analogous interpretation: the input via water from colony percolation, as hypothesised in previous works (Cipro et al., 2018b[START_REF] Cipro | Organic pollutants and their correlation with stable isotopes in vegetation from King George Island, Antarctica[END_REF] would be the more relevant specific exposure phenomenon. This similarity does not occur for soil samples, likely meaning that specific mechanisms at the local scale regulate each one of the situations.

-PCBs chlorination level

The average relative chlorination number of the PBCs found in the matrixes is shown in Figure 1.

In all cases, the profile from the control sites was lighter than the one from colony sites, meaning that lower molecular weight congeners are more subject to undergo a transport process via resuspension/volatilisation. Yet, the results show congeners with 5 and 4 chlorine substitutions as the most frequent, in this order, except for lichen from control sites, when this order was inverted; These results reinforce the hypothesis that the exposure to PCBs only occur through atmospheric deposition due to the lack of a true root system to absorb compounds from the substratum [START_REF] Liu | Response of stable carbon isotope in epilithic mosses to atmospheric nitrogen deposition[END_REF].

Moreover, the comparison with the data from King George Island only, in a study in which the samples were collected further away from any interference of seabird colonies [START_REF] Cipro | Organic pollutants and their correlation with stable isotopes in vegetation from King George Island, Antarctica[END_REF], shows PCBs profiles in vegetation even lighter than the ones from the present study, reinforcing the fact that heavier congeners are less mobile regarding both primary (see von [START_REF] Von Waldow | Quantifying remoteness from emission sources of persistent organic pollutants on a global scale[END_REF] and secondary sources: as one moves away from colonies, the profiles become proportionally lighter due to a lesser input of heavier congeners and also from a higher input of lighter ones. Because heavier congeners are more environmentally persistent [START_REF] Fuoco | Polychlorobiphenyls in Antarctic matrices[END_REF], this has also to be taken into account when considering this proportionally heavier profiles from colonies. The overall profiles for both vegetation subgroups come in close agreement with previous literature, particularly for snow samples (Cipro et al., 2017b), subject to a much similar input mechanism from the atmosphere.

Soil samples, however, concentrated roughly 90% of their congeners in tetra and penta chlorinated congeners, with the latter representing around 75%. This distribution is corroborated by literature [START_REF] Roosens | Penguin colonies as secondary sources of contamination with persistent organic pollutants[END_REF] where pentachlorinated congeners also prevailed. This distribution is due to the fact that most of the more persistent PCBs sit in the interval from 4 to 6 chlorine substitutions. As the authors state, the profiles were dominated by comparatively heavier congeners indicating penguins as the main source rather than long range atmospheric deposition, in which case a lighter profile would be expected .

-Factors influencing PCBs and HCB concentrations

In order to evaluate the relative importance of the categorical variables (main location, specific location, colony species and matrix species), GLM models were built and model selection was based upon AICc as previously explained. Results are presented in Table 3.

No model building was possible for the soil colony subgroup neither for PCBs nor for HCB due to the number of samples under the method detection limit paired with the other variables being taken into account.

The results come partially in agreement with the first part of the study. In a general way, models with less factors were more present among the most accurate in the present study, with models presenting only one factor being the most accurate in almost half of the cases. Specific

Location and %C were the most present factor in the chosen models, whereas for trace elements Specific Location and Matrix Species seem to play a more important role. This means that for POPs, due to the comparatively higher tendency to volatilise/resuspend in regard to trace elements, the organic matter concentration, Matrix Species seem to have a less important role than the sheer quantity of organic matter, except for lichens. This also comes in agreement with previous works suggesting that a migration into the liquid phase is likely to precede the absorption of POPs by more water-dependent organisms (Cipro et al., 2018b[START_REF] Cipro | Organic pollutants and their correlation with stable isotopes in vegetation from King George Island, Antarctica[END_REF] when compared to lichens that will not interact with the substratum and absorb contaminants directly from the atmosphere [START_REF] Yogui | Polybrominated diphenyl ether flame retardants in lichens and mosses from King George Island, maritime Antarctica[END_REF].

-Control sites versus colony sites

Firstly, Tukey HSD/ANOVA tests were performed in order to understand whether some control sites with apparently abnormal concentrations were really not grouped with their counterparts, confirming some hypothesised local sources other than the colonies. Even with apparent huge disparities among the data subsets, the only case when the homogenous groups after the post-hoc test indicated a separation was for PCBs in lichens from Deception Island.

Indeed, this site presented the highest concentrations for PCBs from control sites in both lichen and mosses, but was the least contaminated control soil dataset. This pattern, however, is not repeated for HCB. These results for vegetation are much higher than those found in previous literature [START_REF] Cabrerizo | Influence of organic matter content and human activities on the occurrence of organic pollutants in antarctic soils, lichens, grass, and mosses[END_REF], which presented levels of PCBs in locations closer to ancient stations (such as our location) higher than those found in the vicinity of more recent ones. This means that ancient contamination and eventual spillages might have not yet been completely revolatilised in a homogenous manner as the authors hypothesise [START_REF] Cabrerizo | Influence of organic matter content and human activities on the occurrence of organic pollutants in antarctic soils, lichens, grass, and mosses[END_REF], at least not in our area of sampling within Deception Island, different from theirs. In other words, if a seabird colony can act as a relevant secondary source, an abandoned station also can. Yet, we hereby hypothesise another factor, which is the local volcanic activity via fumaroles, increasing the temperature in some zones, favouring the revolatilisation of PCBs and thus, the exposure of vegetation to them.

Following, the results for the t-tests of colony sites compared to the closest control site are presented in Table 4.

Lichens presented significant differences in POPs concentrations between control and colony sites in two occasions: for PCBs from a chinstrap penguin Pygoscelis antarcticus colony from Half Moon Island and for HCB from an Antarctic tern Sterna vittata colony from Penguin Island. In both cases there has been also a significant difference in the organic matter content, corroborating the statement that the role of nutrient source is not separated from the role of contaminants source.

Mosses, on the other hand, did not present any significant difference in POPs concentrations. Despite some huge differences found, they were not statistically significant likely due to the relatively high standard deviation in these cases. Moreover, in some cases the control sites did not overcome the method detection limits and the proper comparison was not possible.

Soil samples, however, presented significant differences in one occasion: Adélie penguin Pygoscelis adeliae colony samples from Turret Point, King George Island, in similar conditions to the ones previously discussed for lichens: significant differences in organic matter content and also SIA indicating that the nutrient and contaminants sources are indeed the same one. Again, the same issues are repeated: some samples presented elevated standard deviations suggesting higher levels than the control sites, but not in a statistically significant manner.

Comparing the matrixes among them and also with the first part of the study (Cipro et al., 2018a), some trends seem clear. Firstly, the accumulation of contaminants is correlated to the contents of organic matter and SIA results, but in different ways. In lichens, %N seems to be the common factor for the most contaminated samples, whereas in mosses the response is somewhat mixed but δ 13 C seems to correlate in the cases when differences from control to colony sites were the highest for both POPs and trace elements. Soil, in turn, presented %C clearly as the common factor for the highest differences in both studies. However, the interpretation for these specific cases is rather limited and has to be taken with caution. However, it is not possible to affirm that because one proxy correlates to a significant difference, then a high level for such proxy will imply a significant difference in POPs levels between control and colony sites.

All of that considered, a broader, more robust approach for interpreting both POPs and trace elements datasets is clearly given by AICc.

-Conclusions

Seabird colonies act as secondary sources of PCBs and likely of HCB as well. The role of colonies as secondary HCB source is not demonstrated as clear as for PCBs, particularly for lichens, likely due to two reasons. Given the comparatively higher HCB volatility, the 150m distance adopted in order to qualify a site as control might not be fully adequate for this analysis and could be reviewed in the future. The second reason is that some control sites have again presented somewhat high concentrations likely due to local sources other than the colonies themselves. For the other matrixes, especially soil, the role of such colonies as relevant secondary seems clearly demonstrated. Soil seemed the most adequate matrix to study this phenomenon and further more detailed research is needed on the vegetation due to specific uptake mechanisms that could not be deeply taken into account, particularly for lichens. Finally, future studies on population parameters and dynamics could clarify some of the interspecific differences found.

Table 1 .

 1 Mean ± standard deviation results for POPs concentrations in (ng g -1 ww), SIA (‰) and C and N content (%) separated by location, colony species and matrix species. Colony species marked with an * indicate very sparse colonies.

Figure 1 .

 1 Figure 1. Relative average PCBs chlorination number in the sampled matrixes.

  None of the PBDEs nor HCHs has overcome the MDL in none of the samples. Due to a methodology modification regarding a previous work (direct GC-ECD quantification in[START_REF] Cipro | Organic pollutants and their correlation with stable isotopes in vegetation from King George Island, Antarctica[END_REF] versus ECD quantification after GC-MS confirmation in the present study), the remaining POPs could be detected and quantified in only a small fraction of the samples and unfortunately do not allow further discussion: chlordanes were present in detectable amounts in only 5 samples, DDTs in 4, Drins in 4 and mirex in 5 from a total of 150 samples. Therefore, Table1presents the data for HCB, PCBs, SIA and centesimal composition for C and N.

Acknowledgements

Again, the authors would like to acknowledge the financing organisms in the University of La Rochelle: CPER (Contrat de Plan État-Région) and FASEP (Fonds d'étude et d'aide au secteur privé). C.V.Z. Cipro received scholarships from FAPESP (São Paulo Research Foundation, Brazil. Grant #2014/15989-0 and grant #2015/07209-8). M. V. Petry and R. C. Montone received funding from INCT-APA (Instituto Nacional de Ciência e Tecnologia Antártico de Pesquisas Ambientais) (CNPq Process n°574018/2008-5), FAPERJ (E-26/170.023/2008) and were supported by the SECIRM, MMA and MCTI. The IUF (Institut Universitaire de France) is acknowledged for its support to Paco Bustamante as a senior

member. G. Guillou and B. Lebreton from the "Plateforme Analyses Isotopiques" of LIENSs are also to be thanked for their assistance during SIA analyses. Finally, many thanks are dedicated to everyone involved in the fieldwork.

Table 2. Significant Spearman correlations among the studied variables.

Table 3 Parameters composition and relative weights (wi) for the most accurate model for each of the subsets according to AICc.