J Mainka

G Maranzana

J Dillet

S Didierjean

O Lottin

Effect of oxygen depletion along the air channel of a PEMFC on the Warburg diffusion impedance

The usual expression of the finite Warburg element used for analysing electrochemical impedance spectroscopy (EIS) spectra is obtained assuming that the reaction occurs at the electrode/membrane interface and that the oxygen concentration at the gas channel/gas diffusion layer interface is constant. A simple improvement of this expression consists in taking into account the oxygen concentration depletion along the gas channel. This pseudo-2D approach should be more appropriate for the investigation of mass transfer limitations in fuel cell membraneelectrode assemblies (MEA). Starting from experimental data from the literature, numerical simulations show that conclusions about the mass transfer limiting layer and its main characteristics can be significantly modified, which can contribute to a better understanding of oxygen transport in fuel cells. The results also put forward the existence of a critical value of the air stoichiometry below which, close to the air channel exit, no oxygen can access to the active layer. However, the diffusion impedance model does not take account of time dependent current variation effects on the gas concentration, although experimental (12,13) proofs have been brought, that they influence significantly the shape and the size of the measured impedance spectra.

Introduction

Fuel cell development has made considerable progress in recent years, so that their incorporation in mobile and stationary power applications becomes a realistic aim for the near future. Among the different categories of fuel cells (depending on the fuel and on the type of electrolyte), proton exchange membrane fuel cells (PEMFC) are very promising, because of their moderate operating conditions (temperatures lower than 100°C, gas pressures from 1 to 3 atm (1)). Using hydrogen as a fuel, they do not emit any greenhouse gases.

PEMFC fed with air must be able to operate at high current densities with a satisfying efficiency and reasonable gas consumption. Under these conditions, the slow oxygen diffusion to the cathode reaction sites represents a significant impediment, which has a direct influence on the oxygen reduction reaction (ORR) kinetics. Various works highlighted the limiting role of the gas diffusion layer for oxygen diffusion (2,3) in the porous cathode of a PEMFC. In this paper, we investigate the possibilities offered by Electrochemical Impedance Spectroscopy (EIS) for the identification of the main parameters governing the transport of oxygen at the cathode, including the effective diffusivity eff D , the characteristic thickness , and the porosity  of the diffusion (equivalent) medium.

The first one-dimensional steady-state diffusion model proposed by Springer et al. (4) made it possible to analyse fuel cell's performance dependence on water management. Since, many improvements have been brought to the theoretical models. Most of them focus on charge and mass transfer at the cathode (5)(6)(7). For instance, Eikerling and Kornyshev ( 6) developed a model describing potential and concentration gradients through the active layer in the direction perpendicular to the electrode surface. In parallel, EIS allows to investigate the different processes occurring in PEMFCs (8) and to analyse in situ fuel cell's behaviour (9,10). One of the most common expressions of the diffusion impedance in fuel cell equivalent circuits is that of the finite Warburg element (11), which yields results in good agreement with the experimental data, at least qualitatively. However, it is based on simple assumptions: Fick's diffusion, surface reaction, and a constant oxygen concentration at the interface with the gas channel. One way to improve the usual Warburg approach was proposed by Doyle et al. and Guo et al., who built pseudo-2D models for diffusion of lithium in batteries (31,32) in order to estimate the diffusion coefficient Ds of lithium ions in the solid phase of the electrode. In the present work, we propose to improve the expression of the Warburg diffusion impedance by considering a decreasing oxygen concentration profile along the gas channel as an alternative boundary condition. This does not change radically the expression of the diffusion impedance, but has a significant impact on the size and on the shape of the impedance spectra, and as a result, on the identified values of the parameters characterising the transport of oxygen. Note that oxygen depletion can be considered in terms of DC effects and in terms of AC effects. This work focuses on variations of the time-averaged oxygen concentration along the channel only; it does not describe the effects of oxygen concentration AC variations. The importance of these effects was recently pointed out by Schneider et al. (12,13) who showed experimentally that the AC signal applied for impedance measurement was at the origin of the oxygen concentration oscillations along the air channel, with a significant impact on the shape and size of the spectra. AC effects are more important at low frequency and at low air stoichiometry. They should also be linked to some geometrical parameters accounting for the capacity of the channel to store oxygen (its depth for instance) and they can appear even when the mass transfer resistance in the gas diffusion layer is negligible. Accounting for the oxygen concentration time-variation in the air channel would involve a full reconsideration of the diffusion impedance and to our knowledge, there is no analytical expression corresponding to this case in the literature. Therefore, the objective of this work is to put forward the limits of the 1D approach starting from the consideration of the decrease in the oxygen time-averaged concentration along the air channel.

Theoretical Model

Before accessing the active reaction sites, air and hydrogen have to flow through different media: the bipolar plates ensuring the gas supply via the flow-field channels, the gas diffusion layers (GDL) homogenising the gas propagation over the electrodes, and the catalyst layers consisting of a porous mixture of carbon powder and catalyst particles embedded in the electrolyte. Gases have to flow through the pores and the electrochemical reactions occur on (or in) the solid phase. Mass transfer models should take into account as well as possible the geometrical characteristics of these successive media but actually, the finite Warburg element describes only gas diffusion through a homogeneous medium.

Model development

According to the majority of the authors (5,6,14), the oxygen reduction reaction is the main factor limiting the electrical efficiency of PEMFCs, even in optimal conditions. On top of that, the access of oxygen to the reaction sites remains problematic, due to the (possible) presence of liquid water in the MEA and to the low oxygen partial pressure when the fuel cell is fed with air (1). On the other hand, since the cathode is the place of water production, the ionic resistance of the well-hydrated polymer can be neglected. As a consequence, the main hypotheses of the one-dimensional cathode model, where oxygen transport occurs only by diffusion are:

• The fuel cell is assumed isothermal and isobaric.

• Ohmic drops in the active layers and in the GDLs are neglected.

• Mass transfer of O2 and H2O is assumed to occur only by diffusion in the pores of the active layer and of the GDL. • Water diffusion through the solid phase of the active layer is not considered. In order to obtain a pseudo-2D model (Figure 1) that takes into account the decrease in the oxygen concentration along the gas channel, three more hypotheses are made:

• The mass transfer resistance in the gas channel is neglected.

• Time variations of the oxygen concentration in the channel are neglected.

• The Tafel slope b is assumed to be constant between the inlet and the outlet of the air channel.

The gases flow along the channel length L; the electrode surface is given by L × d. With these assumptions, the model can be considered as pseudo-two dimensional, the first dimension being the direction of oxygen diffusion from the channel to the active layer whereas the second dimension is the direction of the channel, with a decrease in the timeaveraged oxygen concentration.

Usually a Butler-Volmer formalism is used to describe oxygen reduction kinetics (15) but for sufficiently high current densities and in steady state, the activation overpotential follows a Tafel law of the form:

 = = [1]
In equation [1], the faradaic current density jf (y) and the oxygen concentration cO2 (x,y) vary along the gas channel, i.e. as functions of the y coordinate (Figure 1). The cell-averaged current density J is given by: ( )

0 1 L f I J j y dy Ld L ==  [2]
By convention, the activation overpotential and the faradaic current density at the cathode side are usually negative. For convenience, we consider their absolute values:

( ),0 f act jy   .
In the following, As mentioned above, mass transfer along the x direction is assumed to occur only by diffusion, which means that the global molar flux density in the x direction

2 2 O O H N N +
has to be null (Figure 1). This is a particular case that happens only when half of the water produced at the cathode is evacuated toward the air channel ( ( )

F y j N f O H 4 / 2 =
) to compensate for the oxygen flux 2 O N in the opposite direction. The other half of the water produced by the cell, as well as the molecules flowing from the anode to the cathode under the effect of the electro-osmotic drag must diffuse through the membrane under the effect of a concentration gradient. Hence, the water flux along the gas channel Note that the question of the validity of the Warburg diffusion impedance in PEM fuel cells does not seem to have been fully addressed yet, considering that these operating conditions are not common: the net water drag coefficient can differ significantly from 0.5 [START_REF] Colinart | Experimental study on water transport coefficient in Proton Exchange Membrane Fuel Cell[END_REF]. The fuel cell being fed by air, the oxygen molar ratio along the GDL/gas channel interface is given by: The inlet concentration of oxygen in dry air being 0 According to the previous hypotheses, mass transport of oxygen through the pores of the diffusion media in a PEMFC is described by Stefan-Maxwell equations (1,17,18). However, since the binary diffusion coefficients of O2/H2O and O2/N2 are close to each other, it seems reasonable to use Fick's 1 st law. Assuming that the electrochemical reaction takes place only at the active layer/membrane interface (x = 0), the 1 st Fick's law can be written:

( ) ( ) ( ) ( ) ( ) ( ) ( ) 02
( ) 2 2 0 O O eff x y c N y D x =  =-  . [8]
The diffusion media being porous, it is necessary to use an effective diffusion coefficient eff D (19,20) taking into account their porosity ε and possibly, the presence of Knudsen diffusion (1,20). The most common expressions obey Archie's law (21,[START_REF] Torquato | Random Heterogeneous Materials: Microstructure and Macroscopic Properties[END_REF]:

m eff D D  = [9]
Where m is an exponent varying between 1.5 and 4. In the case of a 3D medium, there is a wide consensus for using m = 3/2, which corresponds to the differential effective medium approximation introduced by Bruggeman [START_REF] Torquato | Random Heterogeneous Materials: Microstructure and Macroscopic Properties[END_REF][START_REF] Bruggeman | Berechnung verschiedener Physikalischer Konstanten von heterogenen Substanzen[END_REF]. m = 3/2 is probably appropriate for active layers, in which the orientation of the solid phase does not follow a privileged direction (although this is not true in the agglomerate models using a cylindrical description (17,20)). In the gas diffusion layers however, the solid phase can be considered as two-dimensional since the carbon fibres are mainly parallel to the electrodes. In this case, the result of the differential effective medium theory is m = 2 [START_REF] Torquato | Random Heterogeneous Materials: Microstructure and Macroscopic Properties[END_REF]. In a liquid, the diffusion is governed by shocks between the particles and D corresponds to the molecular diffusion coefficient Dmol. In a gas, two cases have to be distinguished depending on the values of the pore size and of the molecule's mean free path  (17):

2 A PN 2 RT    = [10]
Where NA is the Avogadro number and  the mean molecule diameter. If the mean free path is much smaller than the pore size ( << dpore) inter-molecular collisions govern the diffusion process and D = Dmol. If the mean free path is of the same range as the pore size ( ≈ dpore), molecule collisions with the pore wall become significant and the Knudsen diffusion must be considered. This is not the case in the GDL where the typical pore size is between 20 and 50 m (24). On the other hand, Kong et al. (25) showed that in most of the active layers, the pore size is between 30 and 60 nm. Therefore, the Knudsen diffusion has to be taken into account. There exist various more or less sophisticated expressions combining D K and Dmol for estimating the diffusion coefficient D. One of the most simple is the Bosanquet formula

K mol D 1 D 1 D 1 + = .
The Knudsen diffusion coefficient of species i can be expressed by

i pore K i M RT 8 d 3 1 D  = .

DC solution

As the different species are neither consumed nor produced in the diffusion layers, their molar fluxes depend only on the y position. Nitrogen being inert, its flux density is null while the oxygen and water flux densities can be expressed as functions of jf (y):

2 ( ) 0 N Ny = , 2 () () 4 f O jy Ny F =- , ( ) 22 2 f m H O H O jy NN F =+ [11] m O H N 2
stands for the amount of water diffusing through the membrane to the anode; according to the hypotheses leading to equation [4], ( )

2 4 f m HO jy N F =-
. In steady-state, Fick's 2 nd law is given by:

2 02 2 0 eff c D x  =  . [12]
Assuming that the oxygen reduction reaction takes place at the electrode/membrane interface (x = 0), equation [11] can be used as a boundary condition for the oxygen diffusion:

( ) 2 0 4 f O eff x y jy c D xF =  =  . [13]
Knowing the oxygen concentration at the channel/diffusion layer interface ( )

x  = [7] allows to solve the Fick's diffusion equations in steady-state [12,13], which yields the expression of cO2 (x,y):

( ) ( ) ( ) ( ) 02 0 0 02 02 02
, 1 14

y f f eff IS d j y dy jy x c x y c IS H FD   -  = - -  +   . [14]
For comparison, the oxygen concentration in the classical 1D model is expressed by:

( ) 0 02 02 1 4 eff Jx c x c FD    = - -   . [15]
The difference between [14] and [15] lies mainly in the oxygen concentration at the gas channel boundary. In the one-dimensional model, it is constant over the electrode surface, whereas in the two-dimensional model it varies between the gas inlet and the outlet. The pseudo-2D model takes also account of the air stoichiometry SO2 and of the gas humidification via the absolute humidity H at the fuel cell inlet.

AC solution

When a small sinusoidal perturbation Δjf (y,t) is added to the mean (DC) current density ( )

f
f f f f t j y t j y t j y j y iwt  = - =  . [17] ( ) ( ) ( ) ( ) , , ( ) exp
The time-averaged value of the oxygen concentration ( ) Thus, in the expressions of the boundary conditions [7] and [13], the current density is replaced by its time-averaged value ( )

f t
jy. Fick's 2 nd law in AC conditions is given by:

2 02 02 2 eff cc D tx  =  [19]
Replacing the oxygen concentration in the 2 nd Fick's law by [START_REF] Colinart | Experimental study on water transport coefficient in Proton Exchange Membrane Fuel Cell[END_REF], it comes:

( ) ( ) ( ) 02 , cosh sinh eff eff jw jw c x y A y x B y x DD      = +             [20]
Equation [13] yields the boundary condition at the electrode/membrane interface ( 0 x = ), which makes it possible to determine ( ) ( ) ( )

By: ( ) ( ) ( ) ( ) ( ) ( ) ( ) (
2 tanh 4 f eff eff jy jw Ay D F jwD    =-    [22]
Then, cO2(x,y,t) is finally given by:

( ) ( ) ( ) ( ) 2 02 2 2 sinh , , Re , 4 cosh eff f iwt O t eff eff jw x D jy c x y t e c x y F jwD jw D       -              = - +                   [23]
At x = 0, the oxygen concentration variation

    = = -     , [ 24]
Which makes it possible to express the faradaic impedance, defined as:

2 2 0 f act act act O f f f f O f x j c Z j j j c j    =      = = -      . [ 25
]
The first term

act f j   
in [25] is the charge transfer resistance Rct:

( )

() ct f t b Ry jy = , in Ωcm 2 . [26]
The second term in [25] 2 2 0

f act O f O f x j c j c j  =   -   
is the oxygen mass transfer impedance Zco2.

Its expression is obtained from [24] and [17]: The oxygen concentration at x = 0 can be obtained with [18], so that the expression of the diffusion impedance [START_REF] Bird | Transport Phenomena[END_REF] can also be written:

( ) ( ) 2 22 tanh cO d DD d d iw Z y R iw   = , [ 28]
With the low frequency limit Rd 2D :

2 02 4 ( 0, ) D d eff t b R FD c x y  = = , in Ωcm 2 [29]
And the characteristic diffusion time τd:

2 d eff D   =
.

[30]

The expression of the pseudo-2D diffusion impedance is close to that of the onedimensional finite Warburg element (11):

( )

1 tanh d D Wd d iw ZR iw   = , [31] With 1 02 4 ( 0) D d eff t b R FD c x  = = .
[32]

The characteristic diffusion time τd [30] is common to both cases. Actually, the difference lies in the expression of the low frequency limit 1

Comparison of the 1D and pseudo-2D approaches

The pseudo-2D approach gives access to the global and local mass transfer characteristics under specified working conditions. In the following, data from the literature (18,26) are analysed in order to open a discussion about its pertinence and interest.

Local mass transfer

The results that are discussed below come from data published by Bultel et al. (18) for a H2/air PEMFC operating at 60°C with an inlet relative humidity RH of about 90%. The main parameters are summarized in Table I.

Table I:

Values of the parameters used for the analysis of mass transfer (18).

Tcell [°C] 60 RH [-] 0.9 H [-] 0.2 P [atm] 1 SO2 1 [-] 2 (if not mentioned otherwise) A [cm 2 ] 0.5 L 1 [cm] 1 d 1 [cm] 0.5 δGDL [µm] 500 (if not mentioned otherwise) εGDL 0.2 2, 2 O H O GDL eff D [m 2 s -1 ] 3.98 10 -5 x εGDL 2 b [V dec -1 ] 0.22 I [A] 0.25
An explicit expression of the diffusion resistance ( )

2D d
Ry can be easily obtained from [29] and [18]: Ry between the air inlet and the air outlet for air stoichiometries varying from 1 to 8. The curves tend toward infinity for SO2 < 1.44 according to [START_REF] Mainka | Proceedings WHEC[END_REF] and considering the values of the parameters given in Table I.

( ) ( ) 2 2 0 0 2 2 () 4 ( ) 1
There is a critical value of the air stoichiometry ( )

0 02 2 0 02 1 4 crit O eff c S J cH FD  = -+ [34]
In the example of Figure 2, the value of ( )

2 crit O
S  is about 1.44 and for SO2 = 1, there is no oxygen reaching the electrode/membrane interface beyond ylim/L = 0.69. The analysis of the ( )

2D d
Ry profiles makes it possible to determine in which conditions the diffusion impedance can be described by a one-dimensional model (when, for instance, the maximum relative variation over the electrode length

( ) ( ) ( ) ( ) 22 2 2 0 0 DD dd D d D d R L R RL R - = is below 10%).
In the conditions of table I ( 18), a one-dimensional behaviour would be achieved only for air stoichiometries SO2 above 11. Ry between the air inlet and the air outlet for different values of the diffusion layer thickness and for a fixed air stoichiometry SO2 = 2 (starting from data of Bultel et al. (18), Table I).

The curves in Figure 3 confirm that the thickness of the diffusion layer has a significant impact on mass transfer kinetics. When SO2 = 2, reducing its initial value (500 µm (18)) by one half decreases significantly the mass transfer limitations: Another way to estimate the effect of oxygen depletion is to compare the impedance spectra obtained with a usual finite Warburg element [31][32] and those obtained with the expression of the diffusion impedance derived above [START_REF] Fischer | [END_REF][29]. The results in Figure 5 correspond to the operating conditions of table I Figure 5 shows the impedance spectra of a whole fuel cell modelled using a Randles (30) equivalent circuit (Figure 4). This equivalent circuit consists of a membrane Resistance m R , a charge transfer resistance ct R , a double-layer capacity dl C and a diffusion impedance ZW (a finite Warburg element [31] or the diffusion impedance accounting for oxygen depletion along the air channel [START_REF] Fischer | [END_REF][29]33]). The values of the impedances in the equivalent circuit -except the diffusion resistance-are taken from the experimental data of Rubio et al. (26). They are kept constant along the air channel, as well as the characteristic diffusion time The numerical results of Figure 5 are confirmed by the shape of experimental spectra measured along the air channel of a segmented PEMFC (Figure 6). The fuel cell was fed in counter flow with pure hydrogen (SH2 = 1.2) and humidified air (Sair = 3) and operated at a mean current density -2 0.5Acm [START_REF] Mainka | Proceedings WHEC[END_REF]. The experimental data on Figure 6a)

t J = (
show clearly the progressive enlargement of the low frequency loop along the air channel, from the inlet to the outlet. Note that in Figure 6a) and 6b), the low frequency loop cannot be distinguished in the impedance spectra of the firsts segments, probably because of the relatively large high frequency impedance at these locations whereas the (low frequency) mass transfer impedance is the lowest. Starting from expressions [START_REF] Fischer | [END_REF][29], the identification of the local mass transfer parameters [START_REF] Mainka | Proceedings WHEC[END_REF] yields values that are in the typical range of PEMFC's gas diffusion layers. However, the results also show that for quantitatively reliable results, AC effects like those recently put forward by Schneider et al. (12,13) have to be included in the oxygen transport model. Note that in the same way, the possible accumulation of liquid water near the air channel outlet (which should affect also the characteristic diffusion time d  ) is not considered here. Finally, some authors (33) mention the existence of a low frequency inductive loop that can be attributed to Pt oxidation. However, this phenomenon arises usually for frequencies of about 100 µHz or less, which is much lower than the frequency domain investigated in this work.

The simple modifications of the hypotheses leading to the expression of the Warburg impedance show that in usual conditions, mass transfer limitations due to oxygen consumption along the channel can be quite significant: the usual 1D approach remains satisfying only for large values of the air stoichiometry which do not correspond to usual operating conditions.

Mass transport properties of the diffusion medium

EIS can be used to estimate the mass transport properties of the diffusion layers. Since 1 /2 DD d

R

and τd depend on the diffusion medium thickness δ, on the effective diffusivity Deff, and indirectly on the porosity ε [9], these parameters can be identified by fitting experimental impedance spectra. Starting from the expression of the characteristic diffusion time τd [30], of 1 /2 DD d R [32,33] and of the charge transfer resistance ct R [26], one can obtain relations that allow to determine the thickness of the limiting diffusion medium  [35,36]. These expressions depend on the operating conditions and on the impedance parameters τd and 1 /2

DD d R
, which must be identified experimentally: [37]

1 01 02 1 4 D d ct D d JR Fc R    =+   . [35] ( ) ( )
Then, the effective diffusivity eff D and the porosity  corresponding to this value of  can be determined with [30] and [9], respectively. δ and  can be compared to the actual thickness and porosity of the gas diffusion layer and of the active layer, which provides indications about the main origin of oxygen transport limitation. Furthermore, comparing the results of equations [35] and [36] tells whether oxygen transport can be considered as one dimensional or two-dimensional. Examples of results are presented in the following, starting from another set of data measured by Rubio et al. (26) with a single PEM cell of 100 cm 2 operated at 12.5 IA = at 1 atm. The fuel cell was fed by pure H2 and dry air with a stoichiometry SO2 = 2.3 (26). In these conditions, the charge transfer resistance was Rct = 13.2 mΩ and the diffusion resistance was Rd = 7.97 mΩ. Since the cell temperature was not given, it has to be estimated, as well as the air inlet relative humidity. Knowing that the fuel cell was not humidified, we chose Tcell = 60 °C and H = 0.1, which is equivalent to a mean RH equal to 0.5 in the channel. The identified values for δ and Deff are given in Table II. (29,30). Oxygen is assumed to diffuse in gas phase through the GDL and the active layer: the possible presence of liquid is accounted for via the actual value of the porosity ε.

Model δ [m]

Deff [m The differences between the parameters identified using the usual Warburg expression and the one resulting from the pseudo-2D approach show clearly that gas consumption along the air channel has an influence, at least in these operating conditions. The values of the diffusion thickness δ and of the effective diffusion coefficient Deff estimated using the pseudo-2D diffusion impedance are higher: the difference ranges from about 30% for the porosity ε to about 100% for the effective diffusion coefficient. The classical one-dimensional Warburg expression yields a value of the effective diffusion coefficient m 2 s -1 -1.4 10 -6 m 2 s -1 ; Table II]; the corresponding value of the porosity (starting from an equivalent -molecular and Knudsen-diffusion coefficient equal to 6.9 10 -6 × εAL 3/2 [9]) is also within the range generally considered for active layers. However, the estimated diffusion thickness is almost 2 orders of magnitude higher than the typical thickness of an active layer (≈ 10 µm). Hence, a clear conclusion about the limiting layer for mass transfer is not possible with the 1D approach. Using the pseudo-2D approach [37], both the porosity and the equivalent diffusion coefficient are within ranges typical of gas diffusion layers. The estimated thickness

), but this could be the result of gas diffusion in the y direction below the channel rib (which would increase the total diffusion thickness) and/or by a possible mass transfer resistance in the gas channel: the gas diffusion layer seems to be the limiting medium for oxygen diffusion, which is in agreement with some previous works (18).

Conclusions and perspectives

EIS is frequently used to investigate the origin of the mass transfer impedance in cathode gas diffusion and active layers. However, the expression of a finite Warburg element used for analysing impedance spectra is obtained assuming that the oxygen concentration at the gas channel/gas diffusion layer interface is constant: this hypothesis is not valid for usual values of the air stoichiometry, and it can lead to wrong estimates of the mass transfer parameters characterising the diffusion media. An alternative approach is proposed: considering gas consumption along the gas channel entails minor changes to the classical Warburg expression and numerical simulations starting from experimental data from the literature show that conclusions about the mass transfer limiting layer and its main characteristics can be significantly modified. The estimated thickness of the GDL is significantly higher than the actual value, which could be linked to a possible mass transfer resistance in the gas channel and/or to constriction effects below the channel rib.

These first results are mostly of qualitative interest and should be confirmed shortly by specific experiments. On top of that, it should be noted that there still exists a wide range of possible improvements in the analytical expression of the diffusion impedance, among which:

• The consideration of oxygen concentration time-variation in the air channel (12,13). • The consideration of all possible values of the water flux, always assumed equal to ( )

F y j N f O H 4 / 2 =
in the present models.

• The consideration of the 3D structure of the and of the distribution of the oxygen consumption through the electrode.

These issues will be the purpose of future works, which should allow a more accurate estimation of the main mass transfer parameters and of the limiting diffusion layer as functions of the operating conditions and of the position along the air channel. Electrode area (m 2 ) b

Tafel slope (V dec -1 ) cO2

Oxygen concentration (mol m -3 )

Figure 1 .

 1 Figure 1. Schematic representation of the pseudo-2D diffusion model. Oxygen concentration depletion in the channel is taken into account for mass transfer modelling.

2  2 ON 2 O 2 

 2222 denote the fluxes in the y direction (in the air channel, in Mol/s) while denote the flux densities in the x direction (in the gas diffusion layer, in Mol/s/m²). The molar fluxes of oxygen is equal to the molar fluxes of oxygen at the gas channel inlet in O minus the amount consumed along the cathode:

O H 2 

 2 (in the y direction) is a function of y, given by:

  at the air inlet. Assuming that air is an ideal gas and that the mean oxygen concentration at the GDL/channel interface in AC conditions

  signal variations are confined in the diffusion layers implicates that

  the absolute humidity H, which should improve the precision of parameter estimation starting from the global diffusion impedance as well as from its local values along the cell surface.

  Figure 2 depicts profiles of

Figure 2

 2 Figure 2. Profiles of

  function of the diffusion thickness δbelow which the oxygen concentration at the electrode/membrane interface

Figure 3

 3 Figure 3. Profiles of

  by a factor of 1.3 when dividing the diffusion layer thickness by one half.

with an air stoichiometry 2 OS

 2 = 3. The value of the characteristic diffusion timed  = 0.751 s -is taken from an experimental work of Rubio et al. (26), performed in similar conditions. The frequency range is 1 mHz to 10 kHz.

Figure 4 .

 4 Figure 4. MEA equivalent electrical circuit according to Randles (30).

d

  , in order to put forward the effects of variations in the value of cm 2 area single cell. As expected, the impedance spectra do not depend on the mass transfer impedance at frequencies above 10 Hz (5). In the low frequency range however (< 10 Hz), the higher values of the mass transfer resistance ( )2D dRy close to the air channel outlet impact significantly the shape and the size of the spectra. Note that the size of the low frequency loop corresponding to the finite Warburg impedance with 21 e. with dry air) is significantly smaller than that obtained with the local value

Figure 5 .Figure 6 .

 56 Figure 5. Comparison between a PEMFC impedance spectra obtained using a classical finite Warburg element[31] and spectra obtained using the diffusion impedance accounting for oxygen depletion[28-29/33]. The PEMFC is modelled using the Randles equivalent circuit (Figure4). All the electrical parameters, except the diffusion impedance, are taken from experimental data measured by Rubio et al.(26).

  typical for active layers [6.2 10 -7

j

  Exchange current density (A cm -2 ) J Cell current density (A cm -2 ) L Gas channel length (m) Mi Molar weight (kg mol -1 ) NA Avogadro number (mol -1 ) Ni Molar flux (mol m -2 s -1

Table II :

 II Diffusion thickness δ and effective diffusion coefficient identified starting from experimental results of Rubio et al. (26) and comparison with typical values from the literature


Juliette Drupt

Christophe Viel

Claire Dune

Vincent Hugel

ROV localization using ballasted umbilical equipped with IMUs

Keywords: Underwater robotics, Cable model, Localization system

This article describes an affordable and setupfriendly cable-based localization technique for underwater Remotely Operated Vehicles (ROVs), which exploits the piecewise linear shape of the umbilical being equipped with a sliding ballast. Each stretched part of the cable is instrumented with a waterproof Inertial Measurement Unit (IMU) to measure its orientation. Using cable's geometry, the vehicle's location can be calculated in relation to the other fixed or moving end of the cable. Experiments carried out with a robotic system in a water tank prove the reliability of this localization strategy. The study investigates the influence of measurement uncertainties on cable orientation and length, as well as the impact of the IMU location along the cable on localization precision. The accuracy of the localization method is discussed.

I. INTRODUCTION

U NDERWATER communication and localization are particularly challenging since Global Navigation Satellite System (GNSS) radio signals are immediately absorbed by the water. Remotely Operated Vehicles (ROVs) are underwater robots equipped with a cable, namely umbilical, that connects them to the surface to maintain bi-directional real-time communication, supply the robot with energy, but also maintain a lifeline with the robot to avoid losing it.

This communication scheme is essential to accomplish tasks that cannot be fully automated and require precision and adaptability, like wreck exploring or infrastructure inspection and maintenance. However, the umbilical is also subject to entanglement and its shape is difficult to predict, complicating the operator's task. Umbilicals are thus a tradeoff between these advantages and constraints. In line with [START_REF] Viel | Rov localization based on umbilical angle measurement[END_REF], the current work aims at turning the presence of this physical link into an advantage by converting it into a localization system using a model of the umbilical and observing its shape. Georeferenced localization of these robots is particularly important for the deployment and recovery of the robots, but also to locate the observations in the scientific missions.

In the absence of tension, a neutrally buoyant cable takes on an irregular shape, increasing the risk of entanglements with itself or its environment. In order to prevent these risks, a technique mostly used for shallow dives is to hang a ballast at a fixed length on the umbilical to stretch the cable portion between the boat and the ballast. If the ROV is close to the ballast, the cable between them takes the shape of a bell, similar to a catenary curve. However, when the ROV moves away from the ballast, the cable is straightened by the Manuscript submitted in December 2023. tension exerted by the ROV and ballast, and thus conforms into a three-dimensional straight line. A ballast tied to a cable can usually only stretch one part of it, the two parts in some particular configuration. However, a ballast moving freely along the cable, namely sliding, will always find its position at the lowest point, corresponding to its minimum of potential energy, where it stretches both parts of the cable simultaneously. The umbilical can then be assimilated to a piecewise linear model, making it particularly suitable for a cable-based ROV localization strategy.

The contributions of the current work are the following:

• a new cable model that includes the shape of the pulley attached to the ballast,

• an implementation of a ROV's cable-based localization system using both cable-mounted IMUs, and ROVmounted Initial Measurement Units (IMUs) and pressure sensor,

• a study of the influence of measurement errors on model parameter estimation and ROV localization accuracy,

• an experimental validation of ROV localization accuracy in test tanks i) for a fixed attachment point, by varying the position of the IMUs along the cable ii) for a mobile attachment point, attached to a second, manuallycontrolled ROV. Related work on ROV localization and cable-based robot localization is presented in Section II. Section III presents the system under study including assumptions and cable system modelling. Embedded sensor-based model parameter measurement is detailed in Section IV. Experiments are presented in Section V to show the effectiveness of the approach and comments on the results obtained. Discussion and perspectives are exposed in Section VI, before concluding the current work in Section VII.

II. RELATED WORK

The geo-referenced localization of an underwater vehicle is usually obtained by locating the vehicle with respect to a reference point -generally, a surface station -with known GNSS coordinates.

In the last decade, the Underwater Wireless Sensor Networks (UWSNs) have received a lot of attention [START_REF] Awan | Underwater wireless sensor networks: A review of recent issues and challenges[END_REF], [START_REF] Su | A review of underwater localization techniques, algorithms, and challenges[END_REF]. UWSNs consist in several interacting components or nodes including vehicles and sensors which are deployed in a specific area in order to perform collaborative monitoring and data collection tasks. In such networks, the localization of a vehicle node can be deduced with respect to the other components of the system, including anchors with known georeferenced coordinates. An UWSN formed by a set of four anchors geo-localized anchors is studied and implemented in [START_REF] Park | 3d underwater localization scheme using em wave attenuation with a depth sensor[END_REF], [START_REF] Park | 3d underwater localization using em waves attenuation for uuv docking[END_REF]. Each anchor broadcasts ElectroMagnetic Waves (EMW) at specific frequencies, allowing the localization of a an underwater vehicle by using the received signal strength to approximate distances and angles to the anchor nodes. UWSNs are however limited by bandwidth restriction, important propagation delays between the nodes, power constraints, high cost due to the deployment of large sensor sets, as well as the need to calibrate the GNSS position of fixed anchors. In addition, the need to equip an area with fixed anchors restricts sea exploration to places where such sensor sets may be deployed.

Underwater vehicles geo-referenced localization is most usually performed by locating the vehicle with respect to a geo-referenced surface vessel using Ultra-Short Base Line (USBL). USBL positioning systems rely on an array of acoustic transducers placed under a ship and a transponder installed in the underwater vehicle. The pose of the robot can be estimated with respect to the boat by travel time and phase difference. USBL systems are lighter than UWSNs in terms of cost and constraints, since only two engines are implied in the localization process, namely the submarine robot and the surface vessel. USBL are therefore the most widely used systems for underwater vehicles georeferenced localization. These systems are although still limited. Pose estimation accuracy is strongly affected by ocean environmental parameters and by sizing uncertainties of the transceiver array. Such systems are therefore highly dependent on a precise calibration. Pose estimation error increases with the distance between the boat and the vehicle: though USBL systems usually range down to 10,000 m, their ranging error is of 1 or 2% for most common systems, leading to poor vehicle localization accuracy in deep water [START_REF] Liu | Research into the integrated navigation of a deep-sea towed vehicle with usbl/dvl and pressure gauge[END_REF], [START_REF] Zhu | A calibration method of usbl installation error based on attitude determination[END_REF]. These systems are also limited by a low update rate due to travel duration of the acoustic waves. Consequently, USBL systems are generally used in a combination with additional embedded sensors including IMUs, cameras, pressure sensors or Sound Navigation and Ranging (SONAR) [START_REF] Audric | Gaps, a new concept for usbl [global acoustic positioning system for ultra short base line positioning[END_REF], [START_REF] Mandić | Underwater object tracking using sonar and USBL measurements[END_REF].

Alternative approaches rely on vision for localization with respect to a fixed anchor with known GNSS coordinates. Usually, visual methods use characteristic patterns on the anchor for detection and relative localization [START_REF] Beaudoin | Underwater field equipment of a network of landmarks optimized for automatic detection by ai[END_REF], [START_REF] Zhang | U-star: An underwater navigation system based on passive 3d optical identification tags[END_REF]. A limitation of these methods is the need to equip the seafloor with visual beacons, including the calibration of their relative and geo-referenced location. The performance of visual methods is also highly dependent on water turbidity and underwater visibility. These visual methods are therefore less common.

The specificity of the current work is to investigate the localization of ROVsusing the cable that connects it to the surface, without the need for additional equipment on site. If known, the shape of the umbilical cable can determine the position of the ROV at its end point. Several works address the problem of cable shape estimation for tethered robotic system, in terrestrial, aerial and underwater applications. Some cables are designed in order to estimate their shape directly, by the use of proprioceptive sensors. The deformation of optical fiber cables can be measured by reflectometry techniques [START_REF] Duncan | High-accuracy fiber-optic shape sensing[END_REF], [START_REF] Yu | Armless underwater manipulation using a small deployable agent vehicle connected by a smart cable[END_REF], [START_REF] Xu | Review of underwater cable shape detection[END_REF]. Other works consider IMUs placed all along a cable and deduce its shape from IMU measurements integration [START_REF] Frank | Smart tether system for underwater navigation and cable shape measurement[END_REF]. However, these technologies require specific expensive cables or cable instrumentation, which are not the standard in marine robotics. Alternatively, the shape of a cable can be estimated according to a physical or geometrical model. Lumped-massspring [START_REF] Buckham | Dynamics simulation of low tension tethers[END_REF], [START_REF] Gonzã¡lez | Real-time simulation of cable pay-out and reel-in with towed fishing gears[END_REF], [START_REF] Hong | Dynamics modeling and motion simulation of usv/uuv with linked underwater cable[END_REF] or finite-elements [START_REF] Eidsvik | Time domain modeling of rov umbilical using beam equations[END_REF], [START_REF] Blintsov | Development of the mathematical modeling method for dynamics of the flexible tether as an element of the underwater complex[END_REF], [START_REF] Nãžrve Eidsvik | Finite element cable-model for remotely operated vehicles (rovs) by application of beam theory[END_REF] methods can provide a dynamic model of the cable given the forces applied to it, but these methods require a significant amount of computational resources and time and strongly depend on the knowledge of the physical parameters of the system. Alternatively, simpler geometrical models can be used. One of them is the catenary model [START_REF] Laranjeira | Catenary-based visual servoing for tether shape control between underwater vehicles[END_REF], [START_REF] D'antonio | The catenary robot: Design and control of a cable propelled by two quadrotors[END_REF], [START_REF] Drupt | Inertialmeasurement-based catenary shape estimation of underwater cables for tethered robots[END_REF], [START_REF] Drupt | Validity of the catenary model for moving submarine cables with negative buoyancy[END_REF], which gives the shape of a quasi-static, non-neutrally buoyant, hanging cable as a parametric three-dimensional curve. The catenary parameters can be computed from visual observation [START_REF] Laranjeira | Catenary-based visual servoing for tether shape control between underwater vehicles[END_REF] or cable tangent inertial measurements at one or two points [START_REF] Drupt | Inertialmeasurement-based catenary shape estimation of underwater cables for tethered robots[END_REF]. However, this model is only valid if the cable weight or buoyancy is very predominant with respect to the other external and internal forces exerted on the cable [START_REF] Drupt | Inertialmeasurement-based catenary shape estimation of underwater cables for tethered robots[END_REF], [START_REF] Drupt | Validity of the catenary model for moving submarine cables with negative buoyancy[END_REF]. It is therefore not suitable in the presence of currents, for cable with neutral or close-to-neutral buoyancy, or if cable stiffness is too high.

Other tethered systems constrain the cable to a piecewise linear model [START_REF] Mcgarey | Tslam: Tethered simultaneous localization and mapping for mobile robots[END_REF], [START_REF] Viel | Self-management of the umbilical of a rov for underwater exploration[END_REF], [START_REF]Self-management of rov umbilical using sliding buoys and stop[END_REF]. In [START_REF] Mcgarey | Tslam: Tethered simultaneous localization and mapping for mobile robots[END_REF], a piecewise linear shape is deduced from cable angle and tension measurement and entanglements in a terrestrial cluttered environment. In [START_REF] Viel | Self-management of the umbilical of a rov for underwater exploration[END_REF], [START_REF]Self-management of rov umbilical using sliding buoys and stop[END_REF], the cable of a ROV is equipped with a system of ballasts and buoys which stretches the tether to conform it into a simple predictable shape, assimilated to straight lines. In [START_REF] Viel | Rov localization based on umbilical angle measurement[END_REF], the authors present a method to estimate the position of the ROV from this umbilical model using local angles and depth. This method was validated in a pool, using a gold standard measure of the model parameters from an underwater motion capture system that tracked the cable and the ROV. Because underwater motion capture tracking is not available outdoors and in field, the measurement of model parameters from embedded sensors was discussed, and was subject to a preliminary feasibility test, as a proof of concept.

The present work extends the concept exposed in [START_REF] Viel | Rov localization based on umbilical angle measurement[END_REF] to propose a real implementation of the method using embedded sensors to measure parameters, and a more complete model of the umbilical considering the presence of curves at the vicinity of the siding element. The orientation angles of the stretched parts are measured by two watertight IMUs mounted on the cable. The depth of cable ends is given by an embedded pressure sensor. The current approach does not require motorization or a mechanical tether management system, and does not depend on a motion capture system for angle reading. Two cable configurations are studied. In the first configuration, the umbilical connects the ROV directly to a point at the surface. In the second configuration the cable connects a ROV to another one, and corresponds to the use of intermediary ROVs in cable management strategies, namely ROV chains [START_REF] Laranjeira | Catenary-based visual servoing for tether shape control between underwater vehicles[END_REF], [START_REF] Drupt | Inertialmeasurement-based catenary shape estimation of underwater cables for tethered robots[END_REF]. III. SYSTEM MODELING Figure 1 represents a scheme of the system. Cable attachment point on the ROV is denoted E, and the other attachment point is denoted O. The first part of the cable from O is constrained to be vertical down by a fixed anchor placed in A, and a sliding ballast is placed along the cable between A and E. Note that one can take O = A. The sliding ballast is composed of a mass fixed on a pulley, which can slide freeely along the cable. A cartesian reference frame Note that a similar system can be defined by replacing the sliding ballast by a sliding buoy, even though the current work focuses on the ballast configuration. To this end, we introduce the buoyancy parameter s b such that s b = 1 if the sliding element is a ballast and s b = -1 if it is a buoy.

F O = (O, x O , y O , z O ) is

A. Assumptions ans application scope

The following assumptions are considered:

(A1) The ratio of mass to buoyancy of the umbilical is negligible compared to the ballast's weight, or currents applied to the ballast; (A2) The length variation of the umbilical is negligible compared to its length, which is therefore considered constant; (A3) When the umbilical is taut, its geometry can be assimilated to straight lines between defined points, here the ballast and the endpoints of the umbilical; (A4) The umbilical attachment point O and the optional anchor point A are assumed to be strong or heavy enough to be unaffected by the action of the umbilical and the ROV, and so can be considered motionless; (A5) ROV motion and cable length are such that the umbilical remains taut. Constraints to meet this assumption are discussed in [START_REF]Self-management of rov umbilical using sliding buoys and stop[END_REF]. The validity of these assumptions will be discussed further in the experimental part, in Section V. The application scope of the present system is for ROVs having an umbilical shorter than 50 m, in order to ensure assumptions (A1) and (A2). In addition, to respect assumption (A3), the current system requires the umbilical to be flexible and allow the sliding ballast to move freely along it. Therefore, the current method is suitable for the following scenarios:

• Exploration of shallow water from a boat with a depth of less than 50m, with zero or stable currents,

• Ship's hull inspection, navigation under uniform ice, etc., since the sliding ballast must remain submerged without touching any obstacle or the sea floor to keep the cable taut,

• Umbilical between an ROV and its cage in case of deep exploration. The umbilical length is then measured from the cage, which is is considered as an anchor in this study.

The cable between the boat and the cage, then, does not need to meet the assumptions described above.

• Chain of ROVs connected with the same umbilical, the model being applicable for each section of cable linking two consecutive ROVs.

B. Cable model

Let us consider a system composed of a ROV tethered to a point O by an umbilical of length l, as described below. Let us model the pulley of the sliding ballast by a circle C of center C and radius r. Let us define B 1 and B 2 the points where the first and last contact between the cable and the pulley are located. The umbilical is then divided in four parts (Fig. 2):

1) the linear part l 0 = ∥OA∥ between attachment point O and an anchor A 2) the linear part l c1 = ∥AB 1 ∥ between the anchor A and the first contact between the cable and the pulley in B 1 3) the arc made by the pulley l p = ⌢ B 1 B 2 4) the linear part l c2 = ∥B 2 E∥ between the last contact between the cable and the pulley in B 2 and the ROV E.

Let L c denote the length of the cable between the anchor A and the ROV E such that L c = l c1 + l p + l c2 , where L c is fixed. The sliding ballast can move freely on the cable L c . In addition, the cable length l is such that l = l 0 + L c . Note that l 0 = 0 if the anchor A = O, and then one gets l = L c . Let us denote B the intersection of lines AB 1 and B 2 E (see Figure 3). In addition, let us define the following lengths:

l 0 l c1 l c2 l p O A B 1 B 2 E C r C F O
l p1 = ∥B 1 B∥ (1)
l p2 = ∥BB 2 ∥ (2)
l 1 = ∥AB∥ = l c1 + l p1 (3)
l 2 = ∥BE∥ = l c2 + l p2 (4)
L = l 1 + l 2 (5) (6)
and non-oriented angles:

σ = B 1 BC (7)
ψ = B 2 BC (8)
θ = σ + ψ (9)
System parameterization at the vicinity of the pulley is represented in Figure 3.

l p1 l p2 C B 1 B 2 B σ ψ θ C Fig. 3.
Cable shape parameterization at the position of the pulley.

Following the steps described in Appendix A, parameters L, l p1 and l p2 can be expressed from pulley radius r and angle θ = σ + ψ by

l p1 = l p2 (10) = tan π -θ 2 r (11)
L = L c -r (π -θ) + 2 tan π -θ 2 r (12)
Note that lim θ→0 L(θ, r) = +∞, i.e. when cable segments AB 1 and B 2 E become parallel. A saturation can so be added such that if L > L max for a defined L max , one keeps L = L c .

If the pulley is very small, the length L can be simplified to L = L c by neglecting the cable deformation caused by the pulley. The validity and influence of the model of the pulley described in the current section will be discussed in Section V.

Let us define angles α, β, µ and η as follows:

• α is the angle between the z-axis of F O and segment AB around the y-axis of F O , in plane (O, x O , z O ); • β is the angle between segment BE and the z-axis of

F O around the y-axis of F O , in plane (O, x O , z O ); • µ is the angle between segment AB and the z-axis of F O around the x-axis of F O , in plane (O, y O , z O );
• η is the angle between the z-axis of F O and segment BE around the x-axis of F O , in plane (O, y O , z O ).

These angles are represented in Figure 4. In addition, we define l 1x and l 1y as the lengths of the projected vector AB on planes (O, x O , z O ) and (O, y O , z O ) respectively. l 2x and l 2y are defined similarly for vector BE.

l 0 l 1x l 2x O A B E α β F O (a) l 0 l 1y l 2y O A B E µ η F O (b) Fig. 4. Angle parametrization of the system in planes (O, x O , z O ) (4a) and (O, y O , z O ) (4b) Let O E = O x E O y E O z E T denote the coordinates of E in frame F O .
As demonstrated in Appendix B and C, O x E and O y E can be expressed as

O x E = l 1x sin (α) + l 2x sin (β) O y E = l 1y sin (µ) + l 2y sin (η) (13)
where

l 1x = l 1 1 + tan (µ) 2 cos (α) 2 (14)
l 1y = l 1 sin (µ) 2 + cos(µ) cos(α) 2 (15)
l 2x = l 2 1 + tan (η) 2 cos (β) 2 (16)
l 2y = l 2 sin (η) 2 + cos(η) cos(β) 2 (17)
and l 1 and l 2 can be expressed as a function of L, O z E , α, β, µ and η as follows:

l 1 = L -l 2 (18)
l 2 = L cos(α) a1 -s b O z E -l 0 cos(α) a1 + cos(β) a2 (19)
where

a 1 = 1 + tan (µ) 2 cos (α) 2 (20)
and

a 2 = 1 + tan (η) 2 cos (β) 2 (21)
From ( 12), L is a function of θ. As shown in Appendix D, θ can be deduced from angles α, β, µ and η as follows:

θ = arccos cos(α + β) a 1 a 2 - sin(µ) sin(η) a 3 a 4 ( 22
)
where a 1 and a 2 are defined in ( 20) and ( 21), and:

a 3 = sin (µ) 2 + cos (µ) cos (α) 2 (23)
a 4 = sin (η) 2 + cos (η) cos (β) 2 (24)
O x E and O y E can therefore be deduced from variable parameters O z E , α, β, µ and η, and constant parameters L c , r, l 0 and s b , where is s b = 1 in the current case since the sliding element is a ballast. Finally, estimation of the coordinates of E in F O allows estimating the position of the ROV by locating the cable attachment E on the ROV with respect to the other cable attachment point O used as a reference. In the following, E is equally referred to as 'ROV position' and 'cable attachment point on the ROV.

It is worth noticing that, unlike [START_REF] Viel | Self-management of the umbilical of a rov for underwater exploration[END_REF], the current study does not require the knowledge of the forces of the ballast, buoyancy, nor current to estimate the position of the ROV. Indeed, these forces are automatically considered in the models with the orientation and value of the angles.

C. Parameters constraints

The sliding ballast must not enter in collision with an obstacle, the seafloor or the surface, what constrains the minimum depth allowed for the ROV. Let O z floor be the minimum z-coordinate in F O of the seabed or of a submarine obstacle inside the sphere of center A and radius L, which delimits the area in which the ROV can move. First, the anchor must respect the constraint l 0 ≤ O z floor [START_REF] Drupt | Validity of the catenary model for moving submarine cables with negative buoyancy[END_REF] to not touch the seafloor, and so the ballast when it is in contact with the anchor. Then, the ballast will never touch the seafloor if

l 0 + L c ≤ O z floor ( 26
)
This condition can be relaxed if the

z-coordinate O z E of the ROV in F O respects O z E < O z lim ( 27
)
where O z lim depends on the ROV's x and y-coordinates in F O . It has been demonstrated in previous work [START_REF] Viel | Self-management of the umbilical of a rov for underwater exploration[END_REF] that an acceptable value for O z lim is

O z lim = 2 O z floor -h B -L 1 - O x E 2 + O y E 2 L2 ( 28
)
where h B denotes the height of the sliding element. ( 28) can be simplified by:

O z lim = 2 O z floor -h B -L c 1 - O x E 2 + O y E 2 L 2 c ( 29
)
since L c ≈ L with L c < L. The interest of using L c instead of L is that this value is known precisely and is constant.

IV. PARAMETERS MEASUREMENT Section III-B gave the position E of the ROV from constant, calibrated parameters and from the set of variable parameters O z E , α, β, µ and η which need to be measured using embedded sensors. The embedded sensors include an IMU and a water pressure sensor for the ROV and two IMUs positioned along the cable, on segments AB and EB respectively.

A. Measurement of O z E O z E is
O z E = d E -d O (30)
Depth d E is measured using the ROV's water pressure sensor and IMU. The pressure sensor is located in D. The IMU is placed in I r and gives the orientation of the inertial frame F Ir with respect to a world frame F w with its z-axis vertical, upwards. This orientation is represented by a rotation matrix w R Ir . The system is represented in Figure 5. Let us define frame F D with origin D and same orientation as F Ir . The position D E of point E in F D is fixed and is assumed to be calibrated.

E I r F Ir D F D Surface F w S F S
d E = -0 0 1 0 S T D D E (31)
where D E is the position of E in F D in homogeneous coordinates and where S T D is a rigid transformation matrix defined from the rotation matrix w R Ir and the depth at the position of the water pressure sensor d D by the relation:

S T D =     w R Ir 0 0 d D 1     (32)
The depth d D is deduced from the water pressure measurement P D in D by the relation

d D = P D -P 0 ρ water * g (33)
where ρ water is water density of water, g is the value of Earth gravity and P 0 , is the pressure returned by the pressure sensor at zero depth.

If O is a fixed point, d O can be calibrated. Otherwise, for instance if O is placed on a second ROV, d O can be measured similarly to d E .

B. Measurement of angles α, β, µ and η

Angles (α, µ) and (β, η) are measured pairwise by two IMUs placed along the cable. As represented in Figure 6, IMUs indexed 1 and 2 are positioned on segments AB and EB in I 1 and I 2 respectively such that IMU j, j ∈ {1, 2} is tangent to the cable along unit vector Each IMU j measures the orientation of its inertial frame F Ij with respect to the world reference frame F w as introduced in Section IV-A, with z-axis vertical, upwards. The coordinates Ij v ij of each vector v ij in F Ij are fixed and assumed to be calibrated. In addition, let us define the lengths

v ij . O A B E v i1 I 1 l i1 F I1 v i2 I 2 l i2 F I2 F O F w
l i1 = ∥AI 1 ∥ (34) l i2 = ∥EI 2 ∥ (35) ∀j ∈ {1, 2}, the coordinates O v ij can be expressed as O v ij = O R w w R Ij ( 36
)
where the rotation matrix w R Ij is given by the measurements of IMU j and O R w is the rotation matrix between frames F w and F O . Because the choice of the orientation of F O is arbitrary except from it z-axis being vertical, downwards, F O can be selected such that O R w is constant and calibrated, or variable but measured, for instance using an additional IMU.

Finally, angles α and µ can be deduced from O v i1 by the relations

α = angle((I 3 -yy T ) O v i1 , z) (37)
µ = angle((I 3 -xx T ) O v i1 , z) (38)
where I 3 is the 3 × 3 identity matrix and

x = 1 0 0 T (39) y = 0 1 0 T (40) z = 0 0 1 T (41)
Similarly, β and η can be obtained from O v i2 by:

β = angle((I 3 -yy T ) O v i2 , z) (42)
η = angle((I 3 -xx T ) O v i2 , z) (43)
V. EXPERIMENTS Experiments are conducted in pools in order the validate the proposed ROV localization system 1 . Two configurations are studied, including a ROV tethered to a fixed point and a pair of ROVs connected together by a cable (ROV chain). The influence of the position of the IMUs along the cable on ROV localization accuracy is also studied, as well as the influence of model parameters measurement errors and the influence of including the effect of the pulley in the cable model.

A. System

1) Pools: Two different pools have been used for the experiments, denoted indoor pool and outdoor pool respectively. Their dimensions are given in Table I. The cable length used for the experiments was therefore restricted by the dimensions of the pools. The indoor pool is made of steel and the outdoor pool contains big steel parts, what prevents the use of magnetometer data in these pools, leading to IMU yaw drift.

2) Robotic system: The experiments use the BlueROV22 from BlueRobotics, which embeds an IMU, a water pressure sensor and a camera. The cable is equipped with a pair of PhidgetSpatial Precision 3/3/3 High Resolution IMUs3 sealed in waterproof housings. The sliding ballast is composed of a neutrally buoyant pulley holding a 240 g lead mass, with a pulley radius r = 2.175 cm. The sliding ballast and the cable IMUs are pictured in Figure 7.

IMUs measurements are integrated using a Madgwick filter [START_REF] Madgwick | An extended complementary filter for full-body marg orientation estimation[END_REF] in order to compute IMU orientation. In the absence of magnetometer data because of the large steel-made parts in the pools, the IMU yaw drift is approximated as linear and calibrated on a static sequence in order to be corrected from the measurements. The standard deviation on the resulting orientation measurements is evaluated on a static sequence

R = R z (θ y ) R y (θ p ) R x (θ r ) (44)
with yaw θ y , pitch θ p and roll θ r . σ r , σ p and σ y are the standard deviation of the angle measurement, for the roll, pitch and yaw after drift correction respectively.

B. ROV connected to a fixed point

A first experiment is set-up in the indoor pool, featuring a BlueROV tethered to a fixed point by its Fathom Slim communication cable 4 such that l 0 = 0 m and L = 3 m. Cable IMUs are placed such that l i1 = 0.4 m and l i2 = 0.2 m. These positions are selected arbitrarily to be far enough from the cable ends to make sure angle measurement is not disturbed by the cable's stiffness, but also close enough from the ends not to block the displacement of the sliding ballast. A 5-cameras Qualisys underwater motion capture (mocap) system5 is set up in the pool in order to track the three-dimensional system. To this end, the robot, the sliding ballast and the cable are equipped with reflective spherical markers. Cable markers are made by reflective tape stuck around the cable every 20 cm. Tape pieces are thin enough not to affect the movements of the sliding ballast. Markers tracking accuracy is about 1 cm. Motion capture output is used as a reference for evaluating the accuracy of the ROV localization system introduced in the current paper. The system is pictured in Figure 8.

The IMU-based ROV localization described in the current paper is compared with the ground truth ROV location provided by the mocap system, but also to the location computed using mocap measurements of parameters O z E , α, β, µ and η, like in [START_REF] Viel | Rov localization based on umbilical angle measurement[END_REF]. This last estimation is referred to as mocap-based. IMU and mocap measurements are put in the same referential frame F O by a change of coordinate frame, for comparison purpose. The data is processed offline, including a synchronization of clock between the IMU and mocap measurements, as well as the transformation between mocap and sensor frames. Two motion sequences are recorded, indexed 1 and 2. Sequence 1 is a short sequence featuring very slow displacements with small range, while Sequences 2 is a four-minutes sequence showing ROV displacements closer to a real use case, with higher speed, motion variation and displacement range. ROV localization results and error with respect to mocap ground-truth are represented in Figures 9 and10 for Sequences 1 and 2 respectively. The error e is the distance between the ground truth and the estimated ROV position.

One can see that the results of IMU-based ROV localization are particularly close to the ground truth in Sequence 1, with an error of about 0.15 m only. IMU-based ROV localization shows a slightly lower accuracy but the corresponding results are still very close to those of the mocap-based ROV localization, meaning that the cable angles are estimated very accurately by the IMUs. This result shows that the proposed localization method can lead to a very accurate ROV position estimation if the ROV is close to stationary. ROV localization error is higher for Sequences 2. The error is about 0.20 m for the IMU-based approach, which is thus acceptable for the tested configuration. The mocap-based approach results in a 0.02 m error, showing that the accuracy of the method is directly proportional to the accuracy of the angle measurement. These results demonstrate the validity of the theoretical model and the possibility of obtaining a fairly accurate position even when the ROV is in motion, since the umbilical remains stretched by the ballast during the displacement.

C. Influence of IMU positioning

Due to the stiffness of the cable, local cable orientation may vary slightly along the cable and be more or less representative of the global cable segment orientation. The influence of the positions of the IMUs along the cable is therefore examined by reproducing the experiment described in Section V-B with different values of l i1 and l i2 . Two positions are compared for each IMU, with l i1 = 0.40 m or 0.78 m and l i2 = 0.20 m or 0.40 m, resulting in four distinct {l i1 , l i2 } configurations. Sequences 1 and 2 studied in Section V-B feature the {l i1 = 0.40 m, l i2 = 0.20 m} configuration. Three additional sequences indexed from 3 to 5 are recorded in order to test the three remaining {l i1 , l i2 } configurations. Similarly to Sequence 2, Sequences 3, 4 and 5 involve motion close to a real ROV use case, with speed and motion variations and large displacement range.

The resulting ROV localization results and error with respect to mocap ground-truth are represented in Figures 11,[START_REF] Duncan | High-accuracy fiber-optic shape sensing[END_REF] reported in Table III, based on Sequences 2, 3, 4 and 5, where the lowest and highest values for each criterion are highlighed in green and red respectively. Because Sequence 1 is the only quasi-static sequence, it is omitted from this analysis. One can observe that the position of the IMUs appears to have little to no impact on the results. Indeed, all configurations lead to the same error order of magnitude, between 0.15 m and 0.20 m. These results show that the IMU-based method proposed in the current work is therefore flexible with respect to the installation of IMUs on the umbilical, whose position can be chosen so as to obstruct the sliding ballast as little as possible.

-2 -1 0 1 O x R [m] Mocap-based-P IMU-based-P Ground-truth 0 1 2 O y R [m] 0 1 2 O z R [m]

D. ROV chain

A second system is studied, featuring a pair of BlueROVs connected together by a cable equipped with a sliding ballast, and therefore forming a robot chain. The connection cable between the robots is simulated by a 3 m long cord. This system is deployed in the outdoor pool. The main advantage of this pool compared to the indoor pool used in Sections V-B and V-C is that it is larger and deeper, making it easier to deploy a two-robot system. However, it is not compatible with the underwater mocap system used in Sections V-B and V-C, which does not work outdoors. Instead, we take advantage of the embedded cameras of the ROVs to generate a Structure from Motion (SfM) localization baseline. This baseline is assumed to be far less precise than the mocap ground truth used in the previous sequences because of the accuracy of the SfM itself, but also and mostly because transferring the SfM estimated pose to the attachment points involve calibrations which introduce important sources of imprecision. However, an accuracy of a few centimeters is to be expected, making this SfM baseline still suitable for evaluating our cable-based ROV localization method.

-2 -1 0 1 O x R [m] Mocap-based-P IMU-based-P Ground-truth 0 1 2 O y R [m] 0 1 2 O z R [m]
The system is represented in Figure 14, and pictured in Figure 15. ∀ k ∈ {1, 2}, Robot k embeds an IMU in P i r,k , a camera in P c k and a pressure sensor in P d k . Frames F i r,k , F c k and F d k are the frames associated to these sensors, respectively. Frame F O is defined with its x-axis vertical, downwards, and with same yaw angle as F ir,1 . In line with the conclusions of Section V-C, the position of the IMUs along the cable is chosen close to the attachment points, in order to obstruct the sliding ballast as little as possible.

The pose of the camera frames F c1 and F c2 is estimated using the SfM software Colmap [START_REF] Schönberger | Structure-from-motion revisited[END_REF] with respect to a reference frame F wc , which is defined arbitrarily. The position of attachment points O and E in frames F c1 and F c2 respectively is measured at the beginning of the experiment, in order to determine their relative position in F wc according to the SfM reconstruction, along the sequence. Since the transformation between frames F wc and F O is unknown at any time of the sequence, we use the distance between points O and E to characterize the estimation error.

-1 A 60 s sequence is recorded. The distance d between attachment points O and E computed from the SfM reconstruction and from our localization method and the corresponding error are reported in Figure 16, as well as the distance error e d defined as the absolute difference between the reference and estimated distance. One can see that the accuracy is about 0.20 m, which is the same order of magnitude as for the system studied in Sections V-B and V-C, featuring a sliding ballast system tethered to a fixed point. Similarly, we can conclude that the current cable-based ROV localization system can also be used in estimating the relative localization of a pair of ROVs tethered by a cable equipped with a sliding ballast.

0 1 O x R [m] Mocap-based-P IMU-based-P Ground-truth 0 1 2 O y R [m] 0 1 2 O z R [m]

E. Influence of parameters measurement error

Sections V-B, V-C and V-D demonstrate the interest of the proposed cable-based localization method for operational applications, but also that its accuracy is directly impacted by the accuracy of parameters measurements. An analysis of the influence of embedded sensor-based parameter measurement errors and cable length measurement error is therefore provided in the current section.

1) IMU-based angle measurement errors: For any angle γ ∈ {α, β, µ, η}, we define γ(t) the value of the angle at time t and γ(t) the value measured by the embedded sensors for angle θ. The measurement error on θ at time t is then e θ (t) = θ(t)θ(t).

Angle measurement error is computed over Sequences 1 to 5 for IMU-based measurements of angles {α, β, µ, η}, using the ground truth angle value θ reported by the motion capture system. The mean µ θ and standard deviation σ θ for θ ∈ {α, β, µ, η} are given in Table IV. By combining the errors on angles {α, β, µ, η}, the global mean error of angle measurement µ angle and the global standard deviation of angle measurement σ angle are determined. One gets:

µ angle = 0.183 • σ angle = 4.266 •
In order to study the effect of this noise on the ROV position error, we model the IMU-based angle measurement as:

γ(t) = γ(t) + ϵ ( 45
)
where ϵ is a noise. As in [START_REF] Viel | Rov localization based on umbilical angle measurement[END_REF], two noise models are considered: a Gaussian noise such that ϵ = ϵ G where ϵ G ∼ N (µ angle , σ 2 angle ) and a uniform noise such that ϵ = ϵ U where ϵ U ∼ U(µ angle -3σ angle , µ angle + 3σ angle ). The effect of these two noise models on ROV localization is estimated by adding a random noise ϵ G or ϵ U to the ground truth angles measured by the mocap system during Sequences 1 to 5. The mean e and standard deviation σ e obtained for the position error e for each noise model is reported in Table V. These computations use the ground truth value of model parameter O z E . One can see that these values have the same order of magnitude than the error observed with the IMU measurements and reported in Table III. The error corresponding to the uniform noise is approximately twice higher than the one generated with the Gaussian noise, from about 20 cm to 40 cm for the mean value. The error obtained with the Gaussian noise is however closer to the error observed using IMU measurements and given in Table III, with an order of magnitude of about 20 cm for error mean and of 10 cm for the standard deviation. The Gaussian noise models seems therefore more accurate.

O B E v i1 I 1 F I1 v i2 I 2 F I2 F O P ir,1 F ir,1 P ir,2
F ir,2 2) Depth measurement error: As in Section V-E1, the error on depth measurement is characterized in order to study the effect of depth measurement noise on the final ROV position estimation. As for angle measurements, we define O z E (t) the value of O z E at time t, and O ẑE (t) the value measured by the embedded sensors. The error on O z E is denoted e z and defined such that e z (t) = O ẑE (t)-O z E (t). This error is computed on Sequences 1 to 5, leading to a mean value µ z and a standard deviation σ z such that:

P d1 F d1 P d2 F d2 F c1 P c1 F c2 P c2 F wc
µ z = -0.030 m σ z = 0.048 m
As in Section V-E1, the effect of measurement noise on O z E is evaluated using two noise models parameterized with µ z and σ z :

O ẑE (t) = O z E (t) + ϵ G (46)
where

ϵ G ∼ N (µ z , σ 2 z ) and O ẑE (t) = O z E (t) + ϵ U (47)
where the mocap during Sequences 1 to 5 leads to a mean error e and error standard deviation σ e given in Table VI. The angles used in these computations are the ground truth values recorded with the motion capture system. It can be observed that these values are significantly lower than those reported in Table V for angle measurement noise. The effect of depth measurement noise is therefore negligible compared to the effect of IMUbased angle measurement noise. 3) Umbilical length error: The influence of umbilical length measurement error e L is studied, such that the measured length is L = L + e L . Figure 17 represents the impact of an umbilical length measurement error between -0.5 m and 0.5 m, with an ideal measurement of all other parameters which is given by the mocap system. The average error in Sequences 1 to 5 is indicated, as well as error variance. One can observe that the average errors increase proportionally and with the same order of magnitude as the umbilical measurement error: a measurement error of 1 cm will result in a localization error close to 1 cm, which makes this method robust to small measurement error. For the cable of length l = 3 m used in the current experiments, a 1 cm length measurement accuracy can be guaranteed.

ϵ U ∼ U(µ z -3σ z , µ z + 3σ z ).

4) Conclusion on the influence of measurement errors:

In the experimental conditions under study, the IMU-based angle measurement errors have been evaluated to cause a ROV localization error of about 20 cm mean and 10 cm standard deviation for the system under study, with a 3 m-long cable. Depth measurement errors have a lower effect on the errors, with a mean value between 5 cm and 10 cm depending on the noise model used, and about 4 cm standard deviation. Finally, the impact of umbilical length error is significantly lower, leading to only 1 cm error with an measurement accuracy of 1 cm, what is easily achievable for the 3 m-long cable used in the experiments. In conclusion, the main source of errors in the proposed ROV localization system lies in the accuracy of the IMU-based angle measurement. This accuracy may nonetheless be improved in several ways. First, in open water, magnetometer data would be available and allow a more accurate yaw angle estimation. In addition, other IMU measurement filters or IMU fixation systems along the cable may be investigated.

F. Influence of the pulley model

For each of the six experimental sequences considered considered in Section V, we compare the results obtained using the cable model described in the current paper, which includes a model that takes the pulley into account, to the cable model used in [START_REF] Viel | Rov localization based on umbilical angle measurement[END_REF] that does not. These models are denoted P-model and NP-model respectively. The NPmodel can be assimilated to a P-model with a zero pulley radius. Table VII gives the mean, median and standard deviation of the ROV position error e obtained with the P and NP models in Sequences 1 to 5, and Table VIII indicates mean, median and standard deviation of the distance error e d given by the two models on the robot chain sequence. One can see that the P-model is slightly more accurate than the NP-model when using mocap-based parameter measurements, in four sequences out of five. However, the NPmodel slightly outperforms the P-model when using embedded sensor-based parameter measurements. As a result, the Pmodel seems to improve the ROV localization estimation only if the model parameters are known very accurately, but the improvement remains marginal.. When the parameters are measured from the embedded sensor, measurements are not accurate enough to observe any improvement when using the P-model since the main source of imprecision lies in parameter measurement noise, and not in model approximations. While the P-model does not improve the IMU-based localization of the systems considered in Section V, it may however become more valuable for a system featuring a larger pulley, for instance for managing a thicker cable.

VI. DISCUSSION AND PERSPECTIVES

A. Influence of dynamics and waves

The ROV localization method proposed in this paper requires that the cable be kept taut by the ballast, or by the buoy in the reversed configuration. However, each time the robot moves up, down or back, a part of the tether becomes temporarily loose. A similar behavior can be observed in the presence of waves. To counteract these effects and keep the cable stretched, the ballast must be chosen to dive faster than the descent of the ROV or the waves. Since a heavy ballast would induce a significant drag force on the umbilical, it is recommended to limit the robot velocity such as it moves slower than the ballast. These solutions are discussed in more details in [START_REF] Viel | Self-management of the umbilical of a rov for underwater exploration[END_REF].

B. Choice of ballast parameters

To meet Assumption A1, the ballast mass or buoy volume must be chosen to have a faster mechanical response than the cable. A complete study is available in [27, Section 9] and provides a solution according to the umbilical and waves parameters. The drag force induced by the ballast can also be estimated, as described in [START_REF] Viel | Self-management of the umbilical of a rov for underwater exploration[END_REF].

C. Interest of the pulley model

The pulley model introduced in this work gives a more accurate description of the cable shape than previous work [START_REF] Viel | Rov localization based on umbilical angle measurement[END_REF] by taking into account its local deformation at the vicinity of the pulley. However, the use of the pulley model will not improve ROV localization accuracy if measurement error in the model's parameters predominates. In addition, the larger the pulley, the more significant the ROV localization accuracy improvement using the pulley model. The pulley model may also become more worthwhile with more accurate cable angle measurements, which may be obtained, for instance, in open water with magnetometer data available.

D. Geo-referenced localization

If linked to a geo-referenced surface station, the ROV can be geo-referenced using the cable-based localization method presented here, which can be particularly useful for underwater scientific exploration, mapping, or when reaching a place of interest.

E. Limitations of the method

Based on the theory and experiments, some limitations can be observed:

1) As already mention in Section III-A, this method is made for short umbilicals of up to 50 m in length. The main source of ROV localization error lies in the accuracy IMU-based cable angle measurements.

2) The presence of the ballast induces a drag force on the umbilical, and so on the robot. 3) Since the umbilical must be stretched to obtain an accurate localization, strong waves and too fast rotations or displacements of the ROV can disturb the localization. However, some solutions are provided in Section VI-B, and it only takes a few seconds of robot standstill to get a correct estimate again. 4) A collision between the umbilical and an obstacle or the seabed can lead to an erroneous localization by bending the cable so that it no longer conforms to the model. Note that collision with an obstacle of known position and shape can be avoided by the operator since the mathematical model provides the shape of the umbilical. Moreover, a collision can be quickly detected by a discrepancy between the localization and the direction control of the robot. Therefore, this method can be used as an entanglement and collision detector if coupled with another localization system like USBL. 5) The present work only considered short cables of 3 m because of the size constraints of operating in pools. A perspective of this work should be to investigate ROV localization accuracy when using longer cables.

VII. CONCLUSION

This paper considered the localization of a ROV by observing the shape of its umbilical, which is conformed into a predictable shape by a sliding ballast which can move freely along it. This shape is composed of two straight lines connected together by a circular arc that follows the shape of the pulley which holds the ballast. The corresponding cable model is parameterized by the local three-dimensional orientation of the cable at each side of the sliding ballast, which can be measured by two IMUs placed on the cable, and by the relative depth of cable ends.

The cable-based ROV localization strategy was implemented on a real robotic system and evaluated in pools using six sequences, which are representative of real ROV use cases. Five sequences were recorded through a ground truth motion capture system, using a ROV tethered to a fixed point. The sixth sequence was recorded using a SfM-based baseline and a pair of ROVs connected by a cable in the context of a robot chain. The influence of cable-IMU position along the cable was found not to affect cable angles measurement accuracy. The ROV localization accuracy using embedded sensors was reported to be about 20 cm for a 3 m-long cable. The main source of ROV localization error was found to be noise from the IMU-based angle measurement. Finally, the model with the pulley did not improve ROV localization significantly, but it could prove useful with larger pulleys or with greater accuracy of model parameters.

The main assets of the proposed cable-based ROV localization method include its use in open water, its very light, cheep instrumentation and very low computational cost compared to usual acoustic-based techniques. In addition, it provides an estimate of the ROV's position and of the threedimensional shape of its cable. The proposed localization scheme is also independent from the ROV's environment and operation conditions, including stable currents and visual or acoustic disturbances. Future works will investigate the deployment of the cable-based ROV localization system in open water, where more accurate IMU-based cable angle measurements are expected in the absence of magnetic disturbances. These evaluations should be carried out with longer cables to check the influence of cable length on localization accuracy. The current cable-based ROV localization may also be fused with other measurements, in order to provide a more accurate ROV localization estimation. To this end, visual, pressure and inertial measurements may be privileged to preserve the lowcost and low-invasive nature of the ROV localization approach.

Lastly, the proposed localization system may be studied for multiple ROVs tethered together as part of a robot chain-based cable management strategy, where n small additional robots are placed on the cable of a ROV and dedicated to cable shape control (see [START_REF] Laranjeira | Catenary-based visual servoing for tether shape control between underwater vehicles[END_REF], [START_REF] Drupt | Inertialmeasurement-based catenary shape estimation of underwater cables for tethered robots[END_REF]). This cable-based localization estimate may be combined with other sources of positioning information available for only one or a few robots of the chain, for instance a USBL transceiver placed in the ROV at the chain's end.

VIII. ACKNOWLEDGMENT

Same calculation can be made for l 2 2x , l 2 2z using respectively (66)-(68), leading to

s b O z E -l 0 = l 1 cos (α) a 1 - l 2 cos (β) a 2 . ( 75
)
where where a 1 , a 2 , a 3 and a 4 are defined in ( 20), ( 21), ( 23) and ( 24) respectively. θ = ABE can be expressed as:

a 1 = 1 + tan (µ) 2 cos (α)
θ =
where c 2 = -P0 ρwater * g . Depth at E can be deduced from d D and the rotation matrix w R Ir . Let us define frame F S with origin the surface point S located vertically from D and with same orientation as F w (Figure 5).

Depth at E is then:

d E = -S z E (86)
where S z E is the z-coordinate of E in F S . In homogeneous coordinates, one gets:

S D = S T D D D (87)
where S T D is the rigid transformation between F D et F S .

According to the definition of frames F D , F S , F Ir and F w one gets:

S T D =     S R D 0 0 d D 1     =     w R Ir 0 0 d D 1     (88)
and finally:

d E = -0 0 1 0 S T D D E (89)

Fig. 1 .

 1 Figure1represents a scheme of the system. Cable attachment point on the ROV is denoted E, and the other attachment point is denoted O. The first part of the cable from O is constrained to be vertical down by a fixed anchor placed in A, and a sliding ballast is placed along the cable between A and E. Note that one can take O = A. The sliding ballast is composed of a mass fixed on a pulley, which can slide freeely along the cable. A cartesian reference frameF O = (O, x O , y O , z O )is defined with its z-axis vertical, downwards.O

Fig. 2 .

 2 Fig. 2. System model including the pulley of center C. The mass attached to the pulley is not represented.

Fig. 5 .

 5 Fig. 5. Points and frames involved in the computation of d E As demonstrated in Appendix E,

Fig. 6 .

 6 Fig. 6. Points and frames involved in the computation of α, β, µ and η

Fig. 7 .

 7 Fig. 7. Cable IMU in waterproof housing (7a) and 240 g sliding ballast (7b)

Fig. 8 .

 8 Fig. 8. Experimental system. Reflective markers are placed on the robot, the pulley, and along the cable.

Fig. 9 .

 9 Fig. 9. IMU-based, mocap-based and ground truth ROV position and positioning error for Sequence 1 ([l i 1 , l i 2 ] = [0.20,0.40] m). The starting and ending timestamps of the curves differ because the mocap system and sensor recording were not monitored to be simultaneous. Note that the mocapbased and ground truth curves overlap very closely.

Fig. 10 .

 10 Fig. 10. IMU-based, mocap-based and ground truth ROV position and positioning error for Sequence 2 ([l i 1 , l i 2 ] = [0.20,0.40] m). Gaps in the curves computed from the mocap recording are caused by a few system tracking failures. Gaps in the IMU-based estimation are caused by temporary inconsistent measurements from which no localization could be computed.

Fig. 11 .

 11 Fig. 11. IMU-based, mocap-based and ground truth ROV position and positioning error for Sequence 3 ([l i 1 , l i 2 ] = [0.78,0.20] m).

Fig. 12 .

 12 Fig. 12. IMU-based, mocap-based and ground truth ROV position and positioning error for Sequence 4 ([l i 1 , l i 2 ] = [0.78,0.40] m).

Fig. 13 .

 13 Fig. 13. IMU-based, mocap-based and ground truth ROV position and positioning error for Sequence 5 ([l i 1 , l i 2 ] = [0.40,0.40] m).

Fig. 14 .Fig. 15 .

 1415 Fig. 14. Two-agent robot chain with a sliding ballast

Fig. 16 .

 16 Fig. 16. Relative distance estimation result.

Fig. 17 .

 17 Fig. 17. Effect of cable measurement error e L on mean ROV localization error e for e L ∈ [-0.5, -0.2, -0.1, -0.05, -0.03, -0.02, -0.01, 0, 0.01, 0.02, 0.03, 0.05, 0.1, 0.2, 0.5] m.

C. Calculation of l 1 and l 2

 2 The z-coordinate of E in F O is such that:O z E = l 0 + s b l 1x cos (α)s b l 2x cos (β) O z E = l 0 + s b l 1y cos (µ)s b l 2y cos (η)(73)an therefores b O z El 0 = l 1x cos (α)l 2x cos (β)(74)since s b ∈ {-1, 1}. Introducing (69) and (71), one gets

2 a 2 2 Since L = l 1 + l 2 ,Let O x 1 O y 1 O z 1 T

 222121 = 1 + tan (η) 2 cos (β) and O x 2 O y 2 O z 2T denote the coordinates of vectors AB and BE respectively, in frame F O . One gets

  the difference of depth between the cable attachment points O and E, of respective depth d O and d E :

TABLE III STATISTICS

 III ON ROV POSITION ERROR e FOR EACH {l i 1 , l i 2 } CONFIGURATION FOR IMU-based ESTIMATION. THE LOWEST AND

	HIGHEST VALUES FOR EACH CRITERIA ARE HIGHLIGHTED IN GREEN AND
		RED RESPECTIVELY.
	Sequence l i 1 (m) l i 2 (m)	#2 0.40 0.20	#3 0.78 0.20	#4 0.78 0.40	#5 0.40 0.40
	mean (m)	0.175	0.179	0.178	0.191
	median (m)	0.172	0.175	0.159	0.196
	σ (m)	0.057	0.090	0.089	0.064

TABLE VI MEAN

 VI AND STANDARD DEVIATION WITH SIMULATED O z E

	MEASUREMENT ERRORS
	Noise model e (m) σe (m)	ϵ G 0.054 0.034	ϵ U 0.082 0.045

TABLE VII STATISTICS

 VII ON ERROR e FOR P AND NP MODELS FOR SEQUENCES 1-5

	Seq. #1 Seq. #2 Seq. #3 Seq. #4 Seq. #5	mean (m) median (m) σ (m) mean (m) median (m) σ (m) mean (m) median (m) σ (m) mean (m) median (m) σ (m) mean (m) median (m) σ (m)	Mocap-based NP P 0.011 0.011 0.001 0.004 0.004 0.001 0.011 0.007 0.014 0.017 0.015 0.015 0.017 0.014 0.017 0.015 0.012 0.019 0.016 0.015 0.008 0.012 0.011 0.008 0.017 0.013 0.021 0.015 0.008 0.022	IMU-based NP P 0.112 0.110 0.017 0.125 0.124 0.017 0.161 0.156 0.057 0.175 0.172 0.057 0.173 0.168 0.090 0.179 0.175 0.090 0.171 0.150 0.089 0.178 0.159 0.089 0.201 0.206 0.069 0.191 0.196 0.064

TABLE VIII COMPARISON

 VIII BETWEEN THE ERROR e d FOR THE P AND NP MODELS FOR

	THE ROBOT CHAIN SEQUENCE
	mean (m) median (m) σ (m)	IMU-based NP P 0.102 0.093 0.070 0.114 0.110 0.074

  Let P D denote the water pressure measured at D.P D = P water + P athm + P b (82)where P athm is the atmospheric pressure, P b is a sensor bias, and P water is the pressure exerted by the water column, with:P water = ρ water * g * d D (83)ρ water is the volumetric mass of water, g is the value of gravity on Earth and d D is the depth of measurement point D with respect to the surface. Let us define P 0 = P athm + P b , which is the pressure returned by the pressure sensor at zero depth.One can get:P water = ρ water * g * d D -P 0

	arccos( where (•) denotes the Euclidian scalar product. Finally, one AB • BE ) (80) ∥AB∥∥BE∥ gets
		θ = arccos	cos(α + β) a 1 a 2	-	sin(µ) sin(η) a 3 a 4	(81)
	E. Calculation of d E
							(84)
	hence:	d D =	P D -P 0 ρ water * g	=	P D ρ water * g	+ c 2

Video available at https://www.youtube.com/channel/UCi8aAbpu8Hjxw3wBXuQfkHw

https://bluerobotics.com/store/rov/bluerov2/

https://www.phidgets.com/?prodid=32

https://bluerov-solutions.com/produkt/fathom-slim-rov-tether-rov-ready/

https://www.qualisys.com/cameras/underwater/

We acknowledge support from the Laboratoire de COnception de Systèmes Mécaniques et Robotiques (COSMER), France, the Centre National de la Recherche Scientifique (CNRS), France, the Laboratoire des Sciences et Techniques de l'Information, de la Communication et de la Connaissance (Lab-STICC), France and the Ifremer, La Seyne-sur-mer, France.

APPENDIX

A. Calculation of L

As illustrated on the Figure 3, since AB 1 and B 2 E are tangent to the circle C, two rectangular triangles BCB 1 and BCB 2 can be defined, with angles BB 1 B = CB 2 B = π 2 . Let us define angle θ such that:

where

Since BCB 1 and BCB 2 are rectangular, then

and

and so it can be defined

The length of the arc ⌢ B 1 B 2 can be evaluated such that ⌢ B 1 B 2 = l p = r θ, and remark ⌢ B 1 B 2 = 0 when θ = 0, i.e. the cable is horizontal.

One has ∥B 1 C∥ = ∥B 2 C∥ = r because are radius of C. Thus, one gets

Finally, since the rectangular triangles BCB 1 and BCB 2 have a common side BC and ∥B 1 C∥ = ∥B 2 C∥ = r, thus the third sides BB 1 and BB 2 of each triangle have the same length (Pythagore theorem), i.e. ∥BB 1 ∥ = ∥BB 2 ∥. Thus,

Finally, L can be expressed as

Note that when the cable is folded, one has θ = 0 and so L tends to infinity: to solve this particular case, a saturation L > L max for a defined L max is defined. On can take L max = L.

B. Calculation of l 1x , l 2x and l 1y , l 2y

where

In addition,

(67)

Using ( 65) and (67), one gets

In the same way, one can obtain using again (65) and (67) (70)