Simon Benhamou
email: simon.benhamou@cefe.cnrs.fr

Of heading, posture and body rotations derived from data acquired by animal-borne accelerometers, magnetometers and gyrometers, kernel density estimation of the corresponding spherical distributions, and fine-scale movement reconstruction

This paper provides the mathematical expressions, formulated in classical trigonometric terms, that are required (i) to compute an animal's body orientations and rotations in 3D space from data acquired at high-frequency by on-board accelerometers, magnetometers and gyrometers, (ii) to relate 3D rotations to changes in 3D orientation, (iii) to represent their spherical distributions through kernel density estimation (KDE), and (iv) to reconstruct fine-scale movements. The way the huge amount of data that is acquired by these on-board devices can be more easily managed is also considered.

Introduction

Thanks to a number of technological improvements, there has been an increasing interest in deploying animal-borne loggers to investigate how animals move at various scales, and thus in attempting to elucidate the underlying behaviors [START_REF] Williams | Optimising the use of bio-loggers for movement ecology research[END_REF]. In particular, three types of high-frequency (>10 Hz) tri-axial bio-logger are commonly used nowadays to study movements performed at very small scale: accelerometers, magnetometers and gyrometers, which make it possible to assess 3D orientations and rotations.

My aim here is to provide a consistent and rational guide to the treatment of the data acquired by these loggers. A part of this information can be gathered in books on aeronautics or robotics, but is usually expressed in terms of complex mathematical formulae (e.g. involving quaternions). Same or other pieces of information, expressed in terms, more easily understandable by biologists/behavioral ecologists, of trigonometric functions were published in bio-logging oriented papers, but in a rather scattered way and they are not always reliable.

In what follows, I will provide the mathematical expressions, formulated in classical trigonometric terms, that are required (i) to compute an animal's body orientations and rotations in 3D space from data acquired at high-frequency by on-board accelerometers, magnetometers and gyrometers, (ii) to relate 3D rotations to changes in 3D orientation, (iii) to represent their spherical distributions through kernel density estimation (KDE), and (iv) to reconstruct fine-scale movements. I will also consider how to manage the huge amount of data that are acquired by these on-board devices. I will focus mainly on loggers assumed to have been tightly mounted on the animal (usually on its back). Much less information can be obtained when the logger is mounted on a collar, liable to roll around the neck. I will consider it only for reconstructing 2D fine-scale movements of a terrestrial animal.

2

Two frames of reference need to be considered. One is Earth-bound. It is defined by the three orthogonal axes X, Y and Z, which, throughout this paper, will conform the so-called East-North-Up (ENU) convention. The other is animal's body-bounded. It is defined by the surge (U), sway (V) and heave (W) axes, corresponding respectively to the tail→head (posterio-anterior), right→left (transversal), and belly→back (ventro-dorsal) axes (Fig. 1). They are orthogonal to, respectively, the frontal (sway x heave) plane, the sagittal (surge x heave) plane and the coronal (surge x sway) plane (in the few species moving in a stand-up posture, such as humans, the frontal and coronal planes are swapped because the surge axis then corresponds to the back→belly rather than tail→head axis). The direction of the surge axis (tail→head) was chosen to reflect the usual forward movement direction. The directions of the sway (right→left) and heave axes (belly→back) were then defined in this way to have the relative orientations of the three body axes, U, V and W, that match those of the X, Y and Z axes, respectively. The three body directions therefore follow the 'right-hand rule' with the index, middle finger and thumb representing the surge, sway and heave axes respectively. The XY plane, corresponding to the water surface, is horizontal, with the Earth-bound axes X pointing to East and Y pointing to North, while the Earth-bound axis Z points upwards (ENU convention). The animal-bound axes, referred to as the surge, sway and heave axes, correspond to the tail→head, right→left and belly→ back axes, respectively. The 3D orientation is defined by the azimuth, θ, the elevation angle, φ, and the bank angle λ. The azimuth is the horizontal component of heading. It corresponds to the orientation of the surge axis in the XY plane. The elevation angle corresponds to the orientation of the surge axis, measured from the horizontal plane in the vertical plane at azimuth θ. It acts as both the vertical component of heading and the longitudinal component of posture. The bank angle is the transversal component of posture. It is measured from the horizontal level in the animal's frontal plane (i.e. orthogonally to the surge axis) and corresponds to the orientation of the sway axis in this plane. The overall body inclination, ζ, combines the longitudinal (φ) and transversal (λ) components of posture. It corresponds to the deviation of the heave axis from the Z axis. Yaw, α, pitch, β, and roll, ρ, are rotations done around the three main body axes.

Heading corresponds to the surge axis orientation expressed in the Earth-bound frame of reference. In 3D space, it is usually specified by two angular spherical coordinates, arXiv:2310.05820 3 hereafter referred to as the azimuth, θ ∈ [-π, π], measured in the horizontal (XY) plane, and the elevation angle, φ ∈ [-π/2, π/2], measured vertically from the XY plane. Classically, the terms 'azimuth' and 'elevation angle' correspond to the compass bearing and angular elevation above the horizon line of an actual object. By extension, they are used here to specify the heading as the compass bearing and angular elevation of a virtual object towards which the animal would face. The posture corresponds to the body orientation with respect to gravity. It is usually specified in terms of longitudinal and transversal inclinations. The former corresponds to the elevation angle, φ ∈ [-π/2, π/2], which is therefore a component of both heading and posture. The latter, λ ∈ [-π, π], referred to as the bank angle, corresponds to the orientation of the sway axis with respect to the horizontal level in the frontal plane (i.e. orthogonally to the surge axis; Fig. 1). It is sometimes mistaken for the vertical orientation of the sway axis measured from the horizontal (XY) plane, φ' (the way λ and φ' are related is highlighted in §2). Overall, the 3D orientation is therefore unequivocally specified by the three Euler/Tait-Bryan angles (θ, φ, λ). With the ENU convention, one gets θ = 0 when the surge axis points eastwards and θ = π/2 when it points northwards, and φ > 0 and λ > 0 when the surge and sway axes, respectively, point above the horizon line. The azimuth, θ, elevation angle, φ, and bank angle, λ, which correspond to orientations expressed in the Earth-bound frame of reference, or the changes in these orientations, are often weirdly referred to as yaw, pitch and roll, whereas these latter terms are expected to refer to rotations expressed in the body-bound frame of reference (see §6).

In this paper, the key trigonometric relationships will be expressed using rotation matrices. They provide an easy way to determine the new Cartesian coordinates of a vector that was rotated by a given angle ϕ in a fixed frame of reference or, equivalently, that was kept fixed while the frame of reference was rotated by -ϕ. Depending on whether the rotation is done around the first (X or U), second (Y or V) or third (Z or W) axis, these matrices are:

R 1 (ϕ)= 1 0 0 0 cos( ) sin( ) 0 sin( ) cos( )     ϕ - ϕ     ϕ ϕ   , R 2 (ϕ)= cos( ) 0 sin( ) 0 1 0 sin( ) 0 cos( ) ϕ - ϕ         ϕ ϕ   , R 3 (ϕ)= cos( ) sin( ) 0 sin( ) cos( ) 0 0 0 1 ϕ - ϕ     ϕ ϕ      
where, in agreement with the ENU convention, ϕ is measured positively counterclockwise for R 1 and R 3 , but clockwise for R 2 when the rotation axis considered points towards the observer. These matrices are first involved in determining the vector B = (B U , B V , B W ) that applies to the animal's body from the vector D = (D U , D V , D W ) recorded by the device (corresponding to acceleration, geomagnetic field or angular velocity, depending on the type of device considered). The device will be assumed to be tightly mounted on the animal's back in such a way that its surge axis lies in the body's sagittal plane (otherwise, the procedure required to infer B from D is quite challenging; see [START_REF] Johnson | A digital acoustic recording tag for measuring the response of wild marine mammals to sound[END_REF]. However, as the animal's back may not be flat, the device may have been mounted somewhat tilted longitudinally (around the sway axis) by φ * , and/or transversally (around the surge axis) by λ * , leading to

B T = R 2 (φ * )R 1 (λ * )D T , i.e.
arXiv:2310.05820

4

B U = D U cos(φ * ) -(D V sin(λ * )+D W cos(λ * )) sin(φ * ) (1a) B V = D V cos(λ * ) -D W sin(λ * ) (1b) B W = D U sin(φ * ) + (D V sin(λ * )+D W cos(λ * )) cos(φ * ) (1c)
with φ * = sin -1 (-a * U /||a * ||) and λ * = atan 2 (-a * V, -a * W ) [1] , where a * = (a * U , a * V , a * W ) is the acceleration vector recorded by the device mounted on the animal kept level (φ = λ = 0; see §2) before release. Rotation matrices will also apply to the unit vectors that specify orientations in 3D space. These vectors are usually initially expressed in terms of two angular spherical coordinates, ξ ∈ [-π, π] and ψ ∈ [-π/2, π/2]. To be multiplied by a rotation matrix, they have first to be expressed in the equivalent Cartesian form, (cos(ξ)cos(ψ), sin(ξ)cos(ψ), sin(ψ)), which, for simplicity, will be written hereafter with the compact notation [ξ•ψ].

Posture and activity estimation from 3D acceleration data

Accelerometers are sensitive to any type of linear acceleration, including gravity. An animalborne high-frequency tri-axial accelerometer therefore makes it possible to determine both the body posture (i.e. orientation with respect to gravity) and the activity level, once the 'static' (gravity-induced) and 'dynamic' (movement-induced) components have been disentangled.

Call a = (a U , a V , a W ) the acceleration vector recorded at any given time t, where a U , a V , and a W corresponds to the values obtained for the surge, sway and heave axes of the accelerometer, respectively. When this device is at rest, one gets ||a|| = g (9.81 m/s 2 ), while the value obtained for each axis ranges between -g and g, depending on its orientation (-g when the axis considered is kept vertically upwards, g when it is kept vertically downwards, 0 when it is kept horizontal). When the device is borne by a moving animal, the data obtained for each axis corresponds to the combination of a static and a dynamic components. Assuming that dynamic acceleration involves mainly fluctuations that tend to quickly counter-balance each other, the static acceleration at time t can be accounted for by the mean vector a ̅ = (a ̅ U , a ̅ V , a ̅ W ), where a ̅ U , a ̅ V , a ̅ W are the mean values obtained for each axis separately, by averaging the values recorded between times t-∆t/2 and t+∆t/2. The temporal window ∆t. must be long enough to smooth acceleration fluctuations but short enough to keep the best possible resolution. A width ranging between 1 and 3 s (depending on the species; Shepard et al. 2008a) is usually appropriate. The norm of the static acceleration, ||a

̅ || = (a ̅ U 2 +a ̅ V 2 +a ̅ W 2 ) 0.5 ,
should usually remain close to g (otherwise the temporal window ∆t has to be enlarged). As it cannot be mounted at the body mass center, an accelerometer is also subjected to linear accelerations due to body rotations, which may result in underestimating the norm of the static acceleration ||a ̅ || (in contrast, the mean norm a should remain close to g). However, the device is usually mounted sufficiently close to the body mass center to consider that values of a u , a v and a w obtained reflect purely linear movements.

As the device may have been mounted tilted longitudinally by φ * and/or transversally by λ * , the body posture may differ from the device's posture. The static acceleration vector of the [1] atan 2 (y, x) = tan -1 (y/x) ± bπ with b=0 for x>0 and b=1 for x<0 arXiv:2310.05820 5 freely moving animal, A = (A U , A V , A W ), can be derived from the static acceleration vector of the device, a ̅ , based on the relationship A T = R 2 (φ * )R 1 (λ * )a ̅ T , i.e. by applying Eqs (1) with B replaced by A and D replaced by a ̅ . Note that the norm is not affected by the tilt: ||A|| = ||a ̅ |.

The body posture can then be specified by the orientations, measured in the vertical plane from the horizontal plane, of its surge axis, φ = sin -1 (-A U /||A||), sway axis, φ' = sin -1 (-A V /||A||), and heave axis, φ" = sin -1 (-A W /||A||). As these three axes are orthogonal to each other, the body posture can be specified unequivocally by only two angles, which correspond to the two angular spherical coordinates of a unit vector whose Cartesian coordinates are sin(φ), sin(φ') and sin(φ"). The body posture can thus be expressed in terms of overall inclination, ζ = π/2-φ" = cos -1 (-A W /||A||), (2a) which measures the angular deviation between the heave axis and the Z (upward vertical) axis, and postural orientation, χ = atan 2 (sin(φ'), sin(φ)) = atan 2 (-A V , -A U ), (2b) which corresponds to the body-related orientation (forwards-backwards and/or rightwardsleftwards) in which the deviation from the Z axis occurs. This handy representation, based on [χ•φ"] = (sin(φ), sin(φ'), sin(φ")), is comparable to a spherical bull's eye spirit level whose equatorial plane stands for the device's coronal plane and bubble stands for the Z axis direction. It is however often more convenient, for subsequent computations, to specify the posture in terms of elevation angle,

φ = sin -1 (-A U /||A||) = tan -1 (-A U /(A 2 V +A 2 W ) 0.5 ), (3a)
and bank angle, which corresponds to the orientation of the sway axis measured in the frontal plane from the horizontal level, λ = atan 2 (sin(φ'), sin(φ")) = atan 2 (-A V , -A W ),

based on [λ•φ] = (sin(φ"), sin(φ'), sin(φ)). The overall body inclination can be easily derived from the elevation and bank angles through the relationship cos(ζ) = cos(φ) cos(λ).

An accelerometer being an inertial measurement device, any movement-induced acceleration is recorded as the result of an apparent force acting in the reverse direction: backward for a forward acceleration, forward for a backward acceleration (i.e. deceleration), centrifugal for the centripetal acceleration involved in turning. The posture estimates obtained using Eqs (2) and (3) are expressed with respect to a reference plane that is orthogonal to this apparent force. As movement-induced accelerations tend to counter-balance each other over ∆t, only gravity should matter (||A|| ≈ g) and the reference plane remains (roughly) horizontal.

However, this plane will be somewhat tilted when an animal undertakes a strong long-lasting forward, backward or centripetal acceleration (e.g., a bird taking off, landing, or gliding along a circle arc, respectively; e.g. [START_REF] Williams | Can accelerometry be used to distinguish between flight types in soaring birds?[END_REF]. Determining the body posture in this context is quite complex. Just to get a flavor, consider an animal that accelerates in the horizontal plane (i.e., orthogonally to gravity), either tangentially or centripetally. In this simple example, one gets ||A|| ≈ (g 2 +h 2 ) 0.5 , with h = ∆v/∆t for a tangential acceleration, or h = v 2 /r for a centripetal acceleration, where v is the linear speed and r is the radius of curvature. arXiv:2310.05820 6

The reference plane is tilted by tan -1 (h/g) with respect to the horizontal plane, either longitudinally for a tangential acceleration or transversally for a centripetal acceleration. The bank angle has then to be corrected: φ = sin -1 (-A U /||A||) -tan -1 (∆v/(g∆t)) in the former case, λ = atan 2 (-A V , -A W ) ± tan -1 (v 2 /(gr)) (-for a left turn, + for a right turn) in the latter case. Thus, for a turn at constant speed in the horizontal plane, a car with stiff suspensions will show no banking (λ = 0) while it is submitted to a strong transversal acceleration (|A V | = v 2 /r), whereas a bicycle has to adopt the bank angle λ = ±tan -1 (v 2 /(gr)) to turn as expected while it is submitted to no noticeable transversal acceleration (A V ≈ 0). Analyses of the flights of a griffon vulture and an osprey (O. Duriez, unpublished data) showed that, unsurprisingly, the 'bicycle way' applies to soaring birds that loop in a thermal to gain altitude.

The data acquired with the accelerometer can also be used to assess activity level as the

mean of the norm, ||d|| = (d U 2 +d V 2 +d W 2 ) 0.5 , of the vector difference, d = (d U , d V , d W ) = a -a ̅ .
The value of ||d|| will be referred hereafter to as the 'Dynamic Body Acceleration' (DBA). The DBA value computed in this way does not depend on the posture. It can therefore be considered a reliable measure of activity, including when the device is mounted on a collar (i.e. when the device's posture may be markedly different from the body posture). It was initially called (somewhat misleadingly, as it is a scalar) 'vectorial DBA' (veDBA) by [START_REF] Qasem | Tri-Axial dynamic acceleration as a proxy for animal energy expenditure; should we be summing values or calculating the vector?[END_REF], to distinguish it from the previous, mathematically inconsistent, way to combine the three d components, |d U |+|d V |+|d W |, called 'overall DBA' (ODBA; R. [START_REF] Wilson | Moving towards acceleration for estimates of activity-specific metabolic rate in freeliving animals: the case of the cormorant[END_REF], which results in unreliable (posture-dependent) activity values. On average, DBA is proportional to oxygen consumption and therefore can act as a suitable proxy for energetic expenses [START_REF] Gleiss | Making overall dynamic body acceleration work: on the theory of acceleration as a proxy for energy expenditure[END_REF][START_REF] Qasem | Tri-Axial dynamic acceleration as a proxy for animal energy expenditure; should we be summing values or calculating the vector?[END_REF]R. Wilson et al. 2020). It can also provide a proxy for the linear speed (see §8).

Static and dynamic acceleration components, a ̅ and d, have initially to be computed at the same high frequency f as the frequency with which raw data, a, is acquired by using a sliding ∆t-width window, with a ̅ being attributed to the center of the window, which encompasses k = round(f∆t) values. Depending on whether k is odd or even, a ̅ has to be computed by averaging the (k-1)/2 or (k-2)/2 preceding values, the central value, and the (k-1)/2 or (k-2)/2 following values. Nevertheless, the amount of data to process in subsequent analyses can be dramatically reduced (divided by k) by keeping, for each time step ∆t, a single a ̅ through subsampling (i.e. as if a ̅ was computed using a jumping rather than sliding window), and a single DBA value through averaging (see §9).

Heading estimation from 3D accelero-magnetic data

Heading in 3D space is defined by both azimuth, θ, and elevation angle, φ. While the latter only requires a high-frequency tri-axial accelerometer to be estimated (Eq. 3a), the former requires a tri-axial high frequency accelero-magnetometer. Indeed, a tri-axial magnetometer, which records the geomagnetic field in 3D space, makes it possible to compute the azimuth with respect to the magnetic North only if its posture is known. For this purpose, it must be arXiv:2310.05820 tightly associated with an accelerometer (on the same chip, sharing the same axes, and being synchronized with a common clock, but not necessarily working at the same frequency).

A magnetometer must be properly calibrated: a shift of the range of the values obtained (offset), specific to each axis, is expected because the device itself generates a magnetic field (hard iron effect). The sensitivity on each axis may also be differently affected by the local environment (soft iron effect). The raw values must therefore be corrected to obtain calibrated values expressed in the same 0-centered range for the three axes. The range width by itself

does not matter. Call (N U , N V , N W ), (E U , E V , E W ), (S U , S V , S W )
, and (W U , W V, W W ), the raw magnetic values measured for the surge (U), sway (V) and heave (W) axes when the surge axis, kept horizontal, points successively towards magnetic North, East, South and West, respectively. These values must be acquired twice, through two horizontal rotations of the surge axis: (i) with the device that rests on its bottom (a U ≈ a V ≈ 0, a W ≈ -g), and (ii) with the device that rests on its left or right flank (a

U ≈ a W ≈ 0, a V ≈ ±g). The 'bottom rotation' should let the W-indexed values unchanged and result in E U ≈ W U ≈ (N U +S U ) /2 and N V ≈ S V ≈ (E V +W V )/2, whereas the 'flank rotation' should let the V-indexed values unchanged and result in N W ≈ S W ≈ (E W +W W )/2
and in U-indexed values alike those obtained with the 'bottom rotation' (otherwise, swap the labels of the axes). The raw magnetic values, r U , r V and r W can then be translated into the calibrated ones:

m U = (2r U -N U -S U )/(N U -S U ) based on the 'bottom rotation' (4a) m V = (2r V -E V -W V )/(E V -W V ) based on the 'bottom rotation' (4b) m W = sgn(a v )(2r W -E W -W W )/(E W -W W ) based on the 'flank rotation' (4c)
Thus, for any axis rotated horizontally, the m value varies from -1 (when the axis considered points southwards) to 1 (when it points northwards). Rescaling it by Bcos(ϕ), where B is the local geomagnetic field intensity and ϕ is the geomagnetic dip angle (inclination), to determine the actual geomagnetic value is unnecessary because the azimuth estimate rests on the ratio of two calibrated values (see below), which is independent of any rescaling factor. If the identities and directions of the three axes are known, a simpler calibration procedure consists in rotating the device at random (in 3D) for a while to determine the minimum and maximum raw values, r min and r max , for each axis separately. As it is not so easy to explore 3D space evenly, one can alternatively perform only two full turns in the horizontal plane, one with the device set level (a U ≈ a V ≈ 0) to determine r min and r max for the surge and sway axes, and the other with the device tilted on its right, left, front or rear flank (a W ≈ 0) to determine r min and r max for the heave axis. Any raw value r can then be translated into the calibrated one m = s(2r-r min -r max )/(r max -r min ), with s = 1 or s = -1 depending on whether r max is obtained when the axis considered points northwards or southward, respectively.

Call m ̅ = (m ̅ U , m ̅ V , m ̅ W ) the mean calibrated magnetic vector obtained by averaging the magnetic values acquired, for each axis, separately over the short time interval ∆t. As the device may have been mounted tilted transversally by λ * (but this does not affect heading) and longitudinally by φ * , the device's azimuth may differ from the body azimuth when the animal moves banked (see Eq. 19). The body-related mean magnetic vector M = (M U , M V , M W ) can be derived from the device-related one m ̅ as

M T = R 2 (φ * )R 1 (λ * )m ̅ T , i.e.

by applying Eqs (1)

with B replaced by M and D replaced by m ̅ . The body azimuth can then be assessed by combining the body-related static acceleration vector A and mean magnetic vector M. Call M = (M U , M V , M W ) the mean magnetic vector that would have been obtained if the body was leveled while keeping its azimuth θ unchanged:

M T = R 2 (φ)R 1 (λ)M T
, where φ and λ are the body elevation and bank angles (Eqs 3). The magnetic North azimuth expressed with respect to the body surge axis is equal to atan 2 (M V , M U ). The body azimuth, measured positively counterclockwise with respect to magnetic East, is therefore

θ = atan 2 (M U , M V ), i.e. θ = atan 2 ((A V 2 +A W 2 )M U -A U (A V M V +A W M W ), (A V M W -A W M V ) ||A||) (5) An alternative formulation is θ = θ d -atan 2 (c 2 , c 1 )
, where θ d is the azimuth of the device (see Eq. 18), and c 1 and c 2 are derived from the relationship (c 1 , c 2 , c 3 )

T = R 2 (φ)R 1 (λ)[0•φ * ] T .
If necessary, the 'true' azimuth value (i.e. based on geographical East) can be obtained by subtracting the geomagnetic declination. However, for fine-scale movement reconstruction (see §8), this magnetic parameter can also be ignored. In the particular case of an animal that heads up or down (φ ≈ ± π/2, A V ≈ A W ≈ 0), the azimuth (and bank angle) becomes irrelevant: its estimates vary randomly due to noise in a and m measures. Complementary information can then be provided by the orientation of the sway axis,

θ' = atan 2 (M V ||A||, -A U M W ).
By the way, note that the body-related mean magnetic vector M can also be used to compute magnetic 'elevation' and 'bank' angles as sin -1 (-M u /||M||) and atan 2 (-M v , -M w ), which are insensitive to the dynamic accelerations that may bias the computation of the standard (gravity-based) elevation and bank angles, φ and λ. However, the geomagnetic field being characterized by a dip angle ϕ, the reference plane (orthogonal to M) against which the magnetic posture is specified is tilted by π/2-|ϕ| with respect to the horizontal plane, making the magnetic posture hard to interpret. The standard (i.e. gravity-based) posture can yet be assessed using M in the particular case of a soaring bird looping in a thermal [START_REF] Williams | Identification of animal movement patterns using tri-axial magnetometry[END_REF][START_REF] Williams | Vultures respond to challenges of near-ground thermal soaring by varying bank angle[END_REF]. At each loop, the arrowhead of M draws a circle around the Z axis with a radius of ||M||cos(ϕ). The sum ΣM of the M vectors sampled evenly over a loop is therefore vertical and deviates from the heave axis by ζ (ϕ<0) or π-ζ b (ϕ>0), with ζ = cos -1 (|ΣM w |/||ΣM||). The elevation angle, φ = sin -1 (-sgn(ϕ)ΣM u /||ΣM||), should be null (a soaring bird is expected to gain altitude only thanks to the vertical advection of warm air) whereas the bank angle, λ = tan -1 (-sgn(ϕ) ΣM v /|ΣM w |), should be equal to ± tan -1 (v 2 /(gr)), where v is the linear speed and r is the radius of curvature (a soaring bird is expected to turn in bending over as a bicycle does).

Rotation estimations from 3D gyroscopic data

A high-frequency tri-axial gyrometer measures the angular velocity components in the coronal, sagittal and frontal planes. They are referred to as yaw, pitch and roll speeds, α ɺ , β ɺ and ρ ɺ , respectively. To be consistent with the ENU convention, α ɺ and ρ ɺ are positive but β ɺ is negative for a rotation that appears anticlockwise when the corresponding rotation axis points towards the observer (the convention used in the firmware of the device may be different). The body-related angular velocity ( , , ) ρ -β α ɺ ɺ ɺ can be derived from the recorded one

d d d ( , , ) ρ -β α ɺ ɺ ɺ based on the relationship ( , , ) ρ -β α ɺ ɺ ɺ T = R 2 (φ * )R 1 (λ * ) d d d ( , , ) ρ -β α ɺ ɺ ɺ T , i.e. by applying Eqs (1) with B U = ρ ɺ , B V = -β ɺ , B W = α ɺ , D U = d ρ ɺ , D V = -d β ɺ
, and D W = d α ɺ . The body-related rotations yaw, pitch and roll for a given time step ∆t = t i -t i-1 can be computed simply as A high-frequency tri-axial gyrometer is the most reliable device for measuring yaw, pitch and roll. Without it, these three body rotations can be assessed indirectly from data acquired by a high-frequency tri-axial accelero-magnetometer (see §6). Reciprocally, changes of azimuth ∆θ for a given time step ∆t can be estimated using a high-frequency tri-axial accelero-gyrometer (see Eq. 13). However, estimating the current azimuth as the sum of changes of azimuth and initial azimuth (

α i = Σ j , i j α ɺ /f, β i = Σ j , i j β ɺ /f and ρ i = Σ j , i j ρ ɺ /f, for j = 1 ... k,
θ i = 1 i j =
Σ ∆θ j + θ 0 ) quickly accumulates random errors. Hence, accelero-gyrometers may be used in place of accelero-magnetometers to assess azimuth only at very short term. A high-frequency accelero-magneto-gyrometer can be useful to cross-check heading estimates, but involves a non-negligible battery drain. It may help infer high-speed movements with great accuracy (A. [START_REF] Wilson | Locomotion dynamics of hunting in wild cheetahs[END_REF], in particular because the tilting of the reference plane when the animal undertook strong long-lasting accelerations prevents suitable estimation of the posture based only on acceleration data (see §2).

Changes of heading and posture

A given heading can be represented as a point at the surface of a unit sphere centered on the animal's current location and whose reference axes (X, Y, Z) match the geographical ones (East, North, Up). The change of heading is done in the plane defined by the two unit vectors,

(θ i-1 , φ i-1
) and (θ i , φ i ), standing for heading at two successive time steps. It can be represented by the orthodromic arc that links the two corresponding points (Fig. 2), i.e. the minor arc of the 'great circle' (i.e. the largest circle that can be drawn on a sphere) passing through the two points considered, which is the shortest possible surface route (geodesic) between them. The change of heading is therefore characterized by the arc size, ω i ∈ [0, π], and initial inclination,

δ i i ∈ [-π, π],
measured positively counterclockwise from the horizontal level in the animal's frontal plane (i.e. tangentially to the unit sphere) at heading (θ i-1 , φ i-1 ). The heading ( i θ ɶ (q), i φ ɶ (q)) at time t i-1 +q∆t, with q∈[0, 1] and with change of heading assumed to be done at the constant speed ω i /∆t for the short time interval

∆t = t i -t i-1 , is obtained from the relationship [ i θ ɶ (q)• i φ ɶ (q)] T = R 3 (θ i-1 )R 2 (φ i-1 )R 1 (δ i i )[qω i •0] T : i θ ɶ (q) = θ i-1 + atan 2 (cos(δ i i ), cos(φ i-1 )/tan(qω i ) -sin(δ i i ) sin(φ i-1 )) (6a) i φ ɶ (q) = sin -1 (sin(qω i ) sin(δ i i ) cos(φ i-1 ) + cos(qω i ) sin(φ i-1
)) (6b) Reciprocally, the size, ω i , and the local inclination, i δ ɶ (q) with q∈[0, 1], of the change of

heading from (θ i-1 , φ i-1 ) = ( i θ ɶ (0), i φ ɶ (0)) to (θ i , φ i ) = ( i θ ɶ (1), i φ ɶ (1)
) can be computed as:

ω i = cos -1 (sin(φ i-1 ) sin(φ i ) + cos(φ i-1 ) cos(φ i ) cos(θ i -θ i-1 )) (7a)
arXiv:2310.05820 10 i δ ɶ (q) = atan 2 (cos(qω i ) sin(φ i ) -cos((1-q)ω i ) sin(φ i-1 ), cos(φ i-1 ) cos(φ i ) sin(θ i -θ i-1 )) (7b) Eqs ( 6) and ( 7) can also be partly derived from the spherical law of cosines. Note that, while the elevation angle i φ ɶ (q) vary from φ i-1 = i φ ɶ (0) to φ i = i φ ɶ (1), and the local reorientation inclination i δ ɶ (q) varies from initial inclination,

δ i i = i δ ɶ (0), to final inclination, δ f i = i δ ɶ (1), the
product cos( i φ ɶ (q))cos( i δ ɶ (q)), which characterizes the inclination φ max of the reorientation plane for the time step ∆t = t i -t i-1 remains constant (cos(φ max ) = cos( i φ ɶ (q))|cos( i δ ɶ (q))|). Modeled as the orthodromic arc linking these two headings, this reorientation can be suitably represented on a unit sphere in a mixed body-bound/Earth-bound frame of reference: the sphere is centered on the animal's location

x i-1 = (x i-1 , y i-1 , z i-1
), but its reference directions are provided by the X axis pointing eastwards, the Y axis pointing northwards and the Z axis pointing upwards (ENU convention). Azimuths, θ i-1 and θ i , are equivalent to longitudes while elevation angles, φ i-1 and φ i , are equivalent to latitudes. The magnitude of the change of heading, i.e. the arc size ω i , is an 'internal' angle whose vertex corresponds to the center of the unit sphere, while the initial and final arc inclinations, δ

η i i = δ i i -λ i-1 at time t i-1 to η f i = δ f i -λ i at time t i , remains constant (η f i = η i i ) during any change of heading without rolling. Assuming that roll, ρ i = η i i -η f i
, is done at constant speed ρ i /∆t between times t i-1 and t i , the bank angle at any intermediate time t i-1 +q∆t with q∈[0, 1] can be interpolated as

i λ ɶ (q) = λ i-1 + i δ ɶ (q) -δ i i + q(η i i -η f i ) (8)
Finally, it is worth noting that Eq. (7a) reflects that the dot product of two unit vectors is simply equal to the cosine of angular discrepancy between them: cos

(ω i ) = [θ i-1 •φ i-1 ].[θ i •φ i ].
Similarly, the cosine of the size, ϖ i , of the change of posture can be computed as cos

(ϖ i ) = [λ i-1 •φ i-1 ].[λ i •φ i ] = [χ i-1 •φ i-1 "].[χ i •φ i "] = A i-1 A i /(||A i-1 ||||A i ||) (9)
arXiv:2310.05820 11

Relating rotations to changes of heading and posture

In 2D space, the posture is irrelevant and the heading is only defined by the azimuth, so that there is a single possible type of rotation, usually referred to as the turning angle, done in the animal's coronal plane. It therefore corresponds to yaw and is simply equal to the change of heading/azimuth. In 3D space, the change of heading is done in a plane that may be quite different from the animal's coronal plane, and therefore may involve both yaw and pitch. Furthermore, even if some roll does not change the heading, it involves a change in the orientation of the coronal plane with respect to the horizontal plane, and therefore affects how the change of heading can be decomposed in terms of yaw and pitch (Fig. 3). Yaw corresponds to a change of azimuth (

α i = ω i = θ i -θ i-1
) for a turn by α i around the heave axis fixed vertically (φ i-1

= φ i = λ i-1 = λ i = β i = ρ i = 0)
. Pitch corresponds to a change of elevation angle (β i = ω i = φ i -φ i-1 for a turn by β i around the sway axis fixed horizontally without going beyond the vertical (λ i-1 = λ i = 0, θ i = θ i-1 , α i = ρ i = 0). Roll corresponds to a change of bank angle (ρ i = λ i -λ i-1 ) for a turn by ρ i around the surge axis set to a fixed heading (

θ i = θ i-1 , φ i = φ i-1 , α i = β i = ω i = 0
). In the general case where the posture at time t i-1 can be any and where a change of heading from time t i-1 to time t i can generate a passive change (i.e. without roll) of the bank angle (see Eq. 8 with

η i i = η f i ), yaw
, pitch and roll are related to changes of heading/posture in a more complex way.

α i •β i ] T = R 1 (δ i -λ)[ω i •0] T .
In this sequential approach, roll is assumed to occur only before yaw and after pitch (its effect is already included in the λ value) so that the arc inclination expressed with respect to the animal's coronal plane remains constant (equal to δ i λ) during the change of heading. Other formulae apply when yaw, pitch and roll are done simultaneously (see text and Fig. 4). In this example, the values of the different angles α = π/3, β = π/6, ω ≈ π/2.80, δ i -λ ≈ π/5.34 (replace π by 180 to obtain values in degrees). Note that δ i and λ are 'surface' angles measured in the animal's frontal plane, i.e. the plane orthogonal to the surge axis, while α, β and ω are 'internal' angles whose vertex correspond to the center of the unit sphere.

If these tree rotations were done sequentially, with pitch done just after yaw (roll being done only before yaw / after pitch), yaw, α i , and pitch, β i , could be considered the two orthogonal components of heading at time t i expressed in the frame of reference defined by the heading/posture at time

t i-1 : [α i •β i ] T = R 1 (-λ i-1 )R 2 (-φ i-1 )R 3 (-θ i-1 )[θ i •φ i ] T .
They also could be considered the two orthogonal components of the orthodromic arc that represents the change of heading between times t i-1 and t i :

[α i •β i ] T = R 1 (η i )[ω i •0] T , i.e. α i = atan 2 (cos(η i ), cot(ω i )) (10a) β i = sin -1 (sin(η i ) sin(ω i )), (10b) where η i = δ i i -λ i-1 = δ f i -λ i
is the value, kept constant (no concomitant roll), of the reorientation inclination with respect to the animal's coronal plane (Fig. 3). It determines how the change of heading is shared between yaw and pitch, from pure yaw (η i = 0 or |η i | = π) to pure pitch (|η i | = π/2). This sequential decomposition of the change of heading proves useful to model 3D RW in terms of yaw and pitch [START_REF] Benhamou | Mean squared displacement and sinuosity of three-dimensional random search movements[END_REF]. It suitably approximates what happens when rotations are done simultaneously if the change of heading is sufficiently small to get sin(ω i ) ≈ ω i , but generates biased results when ω i is large, as may occur when heading is estimated every ∆t rather than at the high recording frequency f. For example, Eqs (10) cannot account for the perfect balance between yaw and pitch (|α i | = |β i |) that is expected when these rotations are done simultaneously while the reorientation plane is at mid-way between the animal's coronal and sagittal planes (i.e. |η i | = π/4 or |η i | = 3π/4).

The simultaneous occurrence of rotational and translational movement components can be mimicked by virtually subdividing the change of heading/posture, from

(θ i-1 , φ i-1 , λ i-1 ) to (θ i , φ i , λ i ), with roll ρ i = η i i -η f i , in k sub-steps.
The change of heading/posture performed for a given time step ∆t = t i -t i-1 can then be modeled as a series of k blocks of three successive small rotations, done in this order: (1) a rotation around the surge axis by ρ i /k; (2) a rotation around the heave axis by α i,j and (3) a rotation around the sway axis by β i,j , for j = 1 ... k. The concomitant change of location from x i-1 to x i can be accounted for by coupling the rotation ρ i /k around the surge axis with a translation along this axis by 2rsin(0.5ω i /k), which corresponds to the length of the chord of an arc of circle with radius r = 0.5||x i -x i-1 ||/sin(ω i /2) and length rω i /k. The k sub-arcs modeling the change of heading have the same size, ω i /k. Without roll, they are all characterized by the same yaw and pitch decomposition, i.e. α i,j = α i /k and β i,j = β i /k for any j. Otherwise, the incremental change by -ρ i /k of the reorientation inclination with respect to the animal's coronal plane, from η i i to η f i , leads each sub-arc to be decomposed in its own way, which can be derived from Eqs (10): α i,j = tan -1 (cos( i η ɶ (q)) tan(ω i /k)) and β i,j = sin -1 (sin( i η ɶ (q)) sin(ω i /k)), with i η ɶ (q) = i δ ɶ (q)-i λ ɶ (q) = η i iqρ i and q = j/k. When k is sufficiently large, these expressions can be rewritten α i,j ≈ cos( i η ɶ (q)) ω i /k and β i,j ≈ sin( i η ɶ (q)) ω i /k (small-angle approximation). Yaw and pitch for the whole arc can therefore be expressed as

α i = Σ j α i,j = cos( ( )) i q η ɶ ω i and β i = Σ j β i,j = sin( ( )) i q η ɶ ω i . If roll is null, i.e. i η ɶ (q)
= η i for any q, one gets α i = cos(η i ) ω i and

β i = sin(η i ) ω i ,
ρ ɺ = (η i i -η f i )/∆t considered constant for the short time step ∆t, i η ɶ (q) is uniformly distributed between i η ɶ (0) = η i i = η i + ρ i /2 and i η ɶ (1) = η f i = η i -ρ i /2, where η i = (η i i +η f i )/2
is the mean value. Yaw and pitch for the time step ∆t = t i -t i-1 can therefore be inferred as:

α i = (sin(η i i ) -sin(η f i )) ω i /ρ i = 2 cos( η i ) sin(ρ i /2) ω i /ρ i (11a) β i = (cos(η f i ) -cos(η i i )) ω i /ρ i = 2 sin( η i ) sin(ρ i /2) ω i /ρ i , ( 11b
) based on the values of ω i , η i i and η f i derived from the heading/postures (θ i-1 , φ i-1 , λ i-1 ) and (θ i , arXiv:2310.05820 13 φ i , λ i ), estimated from the accelero-magnetic measures made between times t i-1 -∆t/2 and t i-1 +∆t/2 and between times t i -∆t/2 and t i +∆t/2, respectively (Fig. 4).

Time t i-1 Time t i

Azimuth θ i-1 (Eq. 5) Azimuth θ i Elev. Angle φ i-1 (Eq. 3a)

Elev. Angle φ i Bank Angle λ i-1 (Eq. 3b) Bank Angle λ i Fig. 4. Inferring body-related rotations from time series of heading/posture. An animal changes its heading (azimuth and elevation angle) and its posture (elevation and bank angles) through rotations -roll, pitch and yaw -done around its main three body axes -surge, sway and heave, respectively. These rotations can be inferred from the time series of successive headings and postures (expressed in the Earth-bound frame of reference specified by a geographical direction and gravity). The initial and final inclinations of the change of heading expressed with respect to the body coronal plane (η i i and η f i ) provides the interface between the body-bound and Earth-bound frames of reference that is required to back convert changes of orientations into rotations. Recall: [ξ•ψ] stands for (cos(ξ) cos(ψ), sin(ξ) cos(ψ), sin(ψ)).

Reciprocally, the 3D reorientation for the time step ∆t = t it i-1 can be assessed from the actual values of yaw, pitch, and roll derived from angular speeds measured by a highfrequency tri-axial gyrometer (see §4). With a time step ∆t sufficiently short to warrant that the angles obtained in this way remain in the range [-π, π], the size of the change of heading and its mean inclination with respect to the animal's coronal plane can be computed as:

ω i = 0.5ρ i (α i 2 +β i 2 ) 0.5 /sin(ρ i /2) ( 12a
)
η i = atan 2 (β i , α i ) (12b) By coupling postures, (φ i-1 , λ i-1
) and (φ i , λ i ), estimated from the acceleration data recorded between times t i-1 -∆t/2 and t i-1 +∆t/2, and between times t i -∆t/2 and t i +∆t/2, respectively, and rotations α i , β i , and ρ i estimated from gyrometric data recorded between times t i-1 and t i , the change of azimuth between times t i-1 and t i , ∆θ i = θ i -θ i-1 , can therefore be computed as

∆θ i = atan 2 (sin(ω i )(cos(φ i-1 )cos(δ i i )+cos(φ i )cos(δ f i ))/2, cos(ω i )-sin(φ i-1 )sin(φ i )) (13) with δ i i = λ i-1 + η i +ρ i /2 and δ f i = λ i + η i -ρ i /2,
and ω i and η i derived from Eqs. ( 12).

CHANGE OF POSTURE Size

ϖ i = cos -1 ([λ i-1 •φ i-1 ].[λ i •φ i ]) BODY-RELATED ROTATIONS Roll ρ i = η i i -η f i Pitch β i = (cos(η f i )-cos(η i i ))ω i /ρ i Yaw α i = (sin(η i i )-sin(η f i ))ω i /ρ i CHANGE OF HEADING Size ω i = cos -1 ([θ i-1 •φ i-1 ].[θ i •φ i ]) Initial inclination δ i i = Final inclination δ f i = atan 2 (sin(φ i )-cos(ω i )sin(φ i-1 ), cos(φ i-1 )cos(φ i )sin(θ i -θ i-1 )) atan 2 (cos(ω i )sin(φ i )-sin(φ i-1 ), cos(φ i-1 )cos(φ i )sin(θ i -θ i-1 )) η i i = δ i i -λ i-1 η f i = δ f i -λ i
arXiv:2310.05820 14

KDE-based computation of spherical distributions

A useful way to assess a distribution rests on kernel density estimation (KDE; [START_REF] Silverman | Density estimation for statistics and data analysis[END_REF]. KDE generates smoothed histograms that do not depend on the arbitrary choice of the grid origin and bin size, which may bias the representation obtained with standard histograms. Furthermore, building a histogram on a sphere involves a tricky tessellation (R. [START_REF] Wilson | A spherical-plot solution to linking acceleration metrics with animal performance, state, behaviour and lifestyle[END_REF]. To apply KDE to locations on the plane, a virtual horizontal grid with a small cell size is overlaid over the plane and the density is estimated at the center of each cell of the grid by setting an isotropic bivariate kernel, usually Gaussian, on each of the n locations. The same principle applies to the locations on a unit sphere such as the ending points of 3D unit vectors . Note however that this makes sense to compute a distribution (whichever the method used) only when the data is generated by a stationary process. Otherwise, the distribution obtained will arbitrarily depend on when the data considered begins and ends to be collected. It is therefore advised to check whether the mean and the variance obtained for the first, second, and third thirds of the time series of the data considered remain similar [START_REF] Benhamou | Of scales and stationarity in animal movements[END_REF]) before attempting to build the distribution.

µ i = (ξ i , ψ i ) with ξ i ∈ [-π, π] (longitude-like) and ψ i ∈ [-π/2, π
For any location x on the unit sphere, the µ-centered bivariate von Mises-Fisher probability density function with concentration parameter κ is

f(x) = κ exp(κ cos(Ω x.µ ))/(4π sinh(κ)) (14a)
where Ω x.µ is the size of the orthodromic arc xµ (i.e. cos(Ω x.µ ) = x.µ). The mean vector length of this distribution is E(cos(Ω x.µ )) = coth(κ) -1/κ [START_REF] Fisher | Dispersion on a sphere[END_REF]; see also [START_REF] Ahmed | Three-dimensional random walk models of individual animal movement and their application to trap counts modelling[END_REF].

Call ||x-µ|| the straight distance between the mean location µ and any given location x on the unit sphere, which therefore corresponds to the chord length of the arc xµ . One gets cos(Ω x.µ ) = 1 -(x-µ) 2 /2. Furthermore, with κ>6, one gets sinh(κ) ≈ cosh(κ) ≈ exp(κ)/2. The von Mises-Fisher distribution with κ>6 can therefore be suitably approximated as:

f(x) ≈ exp(-0.5(x-µ) 2 /σ 2 )/(2πσ 2 ) (14b) with σ 2 = 1/κ ≈ 1-E(cos(Ω x.µ ))
. Thus, a µ-centered bivariate von Mises-Fisher distribution with concentration parameter κ>6 (i.e. σ<π/2.45) and a zero-centered isotropic bivariate Gaussian distribution with variance σ 2 = 1/κ reflect each other through a µ-centered Lambert azimuthal equal-area projection, which maps orthodromic arcs originating at µ as their corresponding chords. Consequently, the area encompassed within a given cumulative isopleth is the same for both distributions. arXiv:2310.05820 15

The standard-deviation of a distribution used as a kernel in KDE is usually referred to as the smoothing parameter (or bandwidth) h. The value required to reliably estimate a spherical distribution is unlikely to be larger than π/2.45 (with h=π/2.45, the 0.95 cumulative isopleth of a von Mises-Fisher distribution is a circle that delimits a cap corresponding to 1/4 of the sphere area). Consequently, although approximate, Eq. (14b) should act as an appropriate kernel for building most (if not all) KDE-based spherical distributions. The density at the center c = (ξ c , ψ c ) of a given cell of the spherical grid can therefore be estimated as

U c = 2 2 2 1 ( ) 1 exp 2 2 i n i nh h =   - -   π    c µ (15)
This expression is identical to the Gaussian kernel-based expression that applies to planar distributions if one agrees that (c-µ) 2 is the squared Euclidean distance between c and µ: (c-

µ) 2 = (x c -x µ ) 2 +(y c -y µ ) 2 on the plane vs. (c-µ) 2 = (x c -x µ ) 2 +(y c -y µ ) 2 +(z c -z µ ) 2
, with x =cos(ξ) cos(ψ), y =sin(ξ)cos(ψ), and z = sin(ψ), on the unit sphere, i.e.

(c-µ i ) 2 = 2(1 -sin(ψ c ) sin(ψ i ) -cos(ψ c ) cos(ψ i ) cos(ξ i -ξ c )) (16)
When KDE applies to locations on a plane, x i = (x i , y i ) with i = 1 ... n, that are independent of each other and come from an isotropic bivariate Gaussian distribution with standard deviation σ, the best trade-off between a level of smoothing sufficiently high to eliminate noise and sufficiently low to keep relevant information is given by h ref = σ/n 1/6 [START_REF] Silverman | Density estimation for statistics and data analysis[END_REF]). When dealing with 3D orientations µ i = (ξ i , ψ i ), the corresponding value is

h ref = κ -0.5 /n 1/6 = (1-||Σ i [ξ i •ψ i ]||/n) 0.5 /n 1/6 (
see example in Fig. 5). However, as this double assumption never holds true, the h ref value only provides an order of magnitude of the value of smoothing parameter required for reliably representing the distribution. Fig. 5. KDE-based distribution of 3D orientations (i.e. locations on a unit sphere). In this example based on computer simulations, n = 3000 locations (dark green dots) were drawn independently of each other from a bivariate von Mises-Fisher distribution centered on location (0, π/6) with concentration parameter κ = 31.5 (corresponding to standard-deviation σ = π/18). The KDE was obtained using von Mises-Fisher kernels centered on each location with a smoothing parameter h = κ -0.5 /n 1/6 ≈ π/68 (replace π by 180 to obtain values in degrees). The density level decreases from the red area to the pale yellow area. The overall area has been cut at the 0.95 cumulative isopleth. In this particular example, any cumulative isopleth is theoretically a circle, but the low density in the peripheral areas makes them highly sensitive to random sampling fluctuations.

The fraction of distribution attributed to a given c-centered cell of the virtual spherical grid is equal to a c U c , where a c is the cell area. If the grid is based on a regular binning of longitudes-like ξ and latitudes-like ψ with angular bin size υ < h, the cell area is variable: a c = 2υsin(υ/2)cos(ψ c ) ≈ υ 2 cos(ψ c ). A grid with a constant cell area a c = υs can be designed through a regular binning of longitudes-like and sines of latitudes-like, with bin sizes υ and s, respectively. To compute the smallest area encompassed within a given cumulative isopleth (such as the 0.95 and 0.50 ones classically used to determine the overall distribution and its core areas), the cells of the virtual grid involved in the computation have to be sorted in the decreasing order of the densities and the corresponding fractions has to be summed up until the cumulative isopleth value is reached. The portion of area one looks for then corresponds to the set of the cells involved by this summation. When cells have all the same area, fractions and densities are proportional to each other. Sorting and summing can then be both performed using a single variable (density or fraction). When the cell area is variable, the procedure is just a bit more complex: sorting densities and summing the corresponding fractions must be made jointly in parallel rather than through a single operation.

Thus, KDE applies to planar and spherical distributions quite similarly. This holds true for both the location-based (i.e. classical) version (LKDE) and the movement-based version (MKDE). MKDE is more suitable than LKDE when the locations are highly serially correlated, as occurred when dealing with 3D orientation. It extends LKDE by setting kernels (with a variable smoothing parameter) on in-between interpolated locations to simulate random bridges linking successive recorded locations [START_REF] Benhamou | Dynamic approach to space and habitat use based on biased random bridges[END_REF]. With this approach the overall distribution can not only be assessed, but also decomposed in terms of spatial distributions of the frequency of visits and of the mean time spent per visit [START_REF] Benhamou | Beyond the utilization distribution: Identifying home range areas that are intensively exploited or repeatedly visited[END_REF]. A program designed to compute these different distributions for 3D orientations, BRB/MKDE_ORI3D (including LKDE as a particular case), can be freely downloaded from www.cefe.cnrs.fr/fr/recherche/bc/dpb/216-simon-benhamou.

Fine-scale movement reconstruction

The data acquired by high-frequency tri-axial accelero-magnetometers, and possibly a speedometer, makes it possible to infer fine-scale movements between successive relocations of animals that cannot be tracked with a very high frequency. This is in particular the case of aquatic animals such as marine mammals, sea turtles, or penguins, which can be GNSSrelocated only when they come to the surface to breath. For non-aquatic animals tracked over long periods, a relatively low relocation frequency may also have been deliberately chosen to spare the battery of the on-board GNSS receptor. Fine-scale movement reconstruction is often referred to as dead-reckoning [START_REF] Wilson | Dead reckoning: a new technique for determining penguin movements at sea[END_REF]), as it is somewhat similar to the way ancient navigators kept track at sea of the compass-bearing of the harbor.

The successive animal's locations are computed recurrently as

x i = (x i , y i , z i ) = x i-1 + v i ∆t,
where v i is the animal's velocity for the i th elementary time step ∆t. The norm of this vector corresponds to the speed, v i , and its orientation corresponds to the heading (θ i , φ i ). Such a recurrent procedure is however prone to accumulate random errors due to the various noises affecting the speed, magnetic and acceleration measures, from which v i is derived. It can therefore work reliably only for relatively short delays (see [START_REF] Gunner | How often should dead-reckoned animal movement paths be corrected for drift?[END_REF]) between ground-truth, GNSS-based relocations used as 'anchors' to compel the reconstructed movement bouts to start and end at these relocations.

The principle for estimating the heading from tri-axial accelero-magnetic data is detailed in §3. Note however that, here, heading (θ i , φ i ) applies to the step x i-1 →x i rather than to the t i time-stamped location x i . It should therefore be assessed from the accelero-magnetic data recorded between times t i-1 and t i rather than between t i -∆t/2 and t i +∆t/2. The forward speed v i can theoretically be estimated by time-integrating the surge acceleration signal, but this signal is usually far too noisy for this purpose. Despite swimming animals move in a smoother way than terrestrial ones, other methods had to be developed to estimate their speed (with various limitations), e.g. a bending paddle (R. [START_REF] Wilson | Prying into the intimate details of animal lives: use of a daily diary on animals[END_REF], a passive propeller [START_REF] Shiomi | Effect of ocean current on the dead-reckoning estimation of 3-D dive paths of emperor penguins[END_REF] or recording the sound level of the water flow against the body (Le [START_REF] Bras | Three-dimensional space use during the bottom phase of southern elephant seal dives[END_REF]. For terrestrial animals, a possibility would be to use the visual flow from the ground, but the amount of energy it requires would prevent this technique from being usable continuously for extended periods. Fortunately, as the reconstructed movement will be stretched or shrunk to start and end at GNSS-based 'anchors' (see below), one only needs a proxy ṽ i that is proportional to the actual speed v i . Without a device that can act as a speedometer, DBA can be considered a speed proxy, with some limitations. Indeed, for numerous species, DBA appears to be roughly proportional to the speed [START_REF] Bidder | The need for speed: testing acceleration for estimating animal travel rates in terrestrial dead-reckoning systems[END_REF]). However, passive gliding in air or water is characterized by the same very low DBA as resting, while a high DBA value can be associated with a null speed when the animal is shaking its body on the spot. A detailed analysis of the acceleration and magnetic signals might help highlight these active phases involving a null speed.

Movement reconstruction involves merging the velocity data with ground-truth (GNSSbased) relocations used as 'anchors' to specify the start and end of each reconstructed movement bout. This procedure removes the systematic errors that may affect the azimuth (e.g. due to neglecting the magnetic declination) and the speed (due to considering only a speed proxy proportional to the actual speed). The accelero-magnetometer (and speedometer, if any) measures are usually time-stamped with an independent clock. If so, these measures must be synchronized with GNSS relocations manually, to the closest second, e.g. by inducing a short and strong acceleration at a precise GNSS time, at the beginning of the tracking period, and if possible (battery still alive), at its end, so as to correct for a possible a time drift of the accelero-magnetometer clock. For marine animals breathing at the surface and equipped with a time-depth recorder (TDR) sharing the same clock as the acceleromagnetometer, the synchronization with GNSS locations can also be done based on surfacing events, as no relocation can be recorded when the animal is underwater. I will first detail the general approach that applies to 3D movements. The simplified version that applies to 2D movements of terrestrial animals for which the logger is mounted on a collar rather than tightly mounted on the animal's back will presented afterwards.

Consider the overall movement defined by a series of N GNSS-based relocations, expressed in terms of longitudes (Ξ I ), latitudes (Ψ I ), and altitudes or depths (Z I ), recorded at times T I , for I = 0 ... N-1. For animals that move over limited extents, longitudes and latitudes can be converted into Cartesian coordinates, (X I , Y I ), through a small-extent metric projection (e.g. UTM). Otherwise, the arc-step linking any two successive relocations (Ξ I-1 , Ψ I-1 ) and (Ξ I , Ψ I ) is usually sufficiently short to be suitably converted into a vector-step through a simple local flat-Earth approximation. In this case, a convenient solution for reconstructing the movement bout between times T I-1 and T I consists in (i) setting X I-1 = 0, Y I-1 = 0, X I = R(Ξ I -Ξ I-1 )cos((Ψ I-1 +Ψ I )/2) and Y I = R(Ψ I -Ψ I-1 ), where R = 6,371,000 m is the mean Earth radius and longitudes and latitudes are expressed in radians, (ii) computing the dead-reckoned Cartesian coordinates x i and y i of in-between locations (see below), and (iii) back-converting This step can also be assessed as the dead-reckoned vector

x i into longitude, ξ i = Ξ I-1 +x i /(Rcos((Ψ I-1 +Ψ I )/2)),
Ṽ I ∆t[ I Θ • I Φ ] = Σ i ṽ i ∆t[θ i •φ i ], for i = 1 ... n = round((T B -T A )/∆t)
, where ṽ i , θ i and φ i are the speed proxy, body azimuth and body elevation angle estimated for the i th elementary time step ∆t:

Ṽ I = (Σ 2 i ṽ xi +Σ 2 i ṽ yi +Σ 2 i ṽ zi ) 0.5 , I Θ =
atan 2 (Σ i ṽ yi , Σ i ṽ xi ) and I Φ = sin -1 (Σ i ṽ zi /Ṽ I ), with ṽ xi = ṽ i cos(φ i ) cos(θ i ), ṽ yi = ṽ i cos(φ i ) sin(θ i ) and ṽ zi = ṽ i sin(φ i ). Systematic errors and the accumulation of random errors result in a deadreckoned vector that does not match the ground-truth vector. As all these errors come from misestimates in terms of orientation and length, the dead-reckoned movement must be corrected in polar (rather than Cartesian, as proposed by R. [START_REF] Wilson | Prying into the intimate details of animal lives: use of a daily diary on animals[END_REF]) terms, through rotation and rescaling (i.e. stretching or shrinking), to preserve its overall shape.

The orientation error can be corrected by rotating the dead-reckoned vector first by -I Θ around the X axis (i.e. in the horizontal plane) to get its azimuth pointing eastwards, then by Φ I -I Φ around the Y axis (i.e. in the E-W vertical plane) to correct its elevation angle without changing its azimuth, and finally by Θ Ι around the X axis (i.e. in the horizontal plane) to fully align it onto the ground truth vector. The error in length can be corrected by rescaling the dead-reckoned vector by L I /(Ṽ I ∆t). The orientation adjustment and rescaling apply similarly to every step. With x 0 set to X I-1 , the series of subsequent locations x i = (x i , y i , z i ), for i = 1 ... n I , specifying the path performed between the two successive ground-truth relocations, can therefore be computed recurrently as

x i T = x i-1 T + R 3 (Θ I )R 2 (Φ I -I Φ )R 3 (-I Θ ) ṽ i T L I /Ṽ I , i.e. x i = x i-1 + (r i cos(Θ I ) -q i sin(Θ I )) L I /Ṽ I (17a) y i = y i-1 + (q i cos(Θ I ) + r i sin(Θ I )) L I /Ṽ I (17b) z i = z i-1 + (ṽ zi cos(Φ I -I Φ ) + p i sin(Φ I -I Φ )) L I /Ṽ I ( 17c
) with p i = ṽ xi cos( I Θ ) + ṽ yi sin( I Θ ), q i = ṽ yi cos( I Θ ) -ṽ xi sin( I Θ ) and r i = p i cos(Φ I -I Φ ) -ṽ zi sin(Φ I -I Φ ). This formulation warrants that the dead-reckoned movement between times T I-1 and T I starts at X I-1 and ends at X I . The location errors at intermediate times tend to increase arXiv:2310.05820 19 with the delay T I -T I-1 , which should therefore be kept as small as possible, in particular when the speed is poorly estimated. For marine animals that breathe at the sea surface (Z I-1 = Z I = 0), Eqs (17) apply with Φ I = 0 and L I = ((∆X I ) 2 +(∆Y I ) 2 ) 0.5 . Furthermore, the dead-reckoned z i values should be replaced by TDR-based ones, if available, as they are more accurate. It is also worth noting that the speed proxy of marine animals applies to the swimming, i.e. water masses-related movement, whereas the GNSS-based relocations used as 'anchors' are expressed in a ground-bound frame of reference. When the water speed is not negligible with respect to the animal's swimming speed, the ground-related vectors linking successive GNSS relocations should first be translated into water-related ones (through vector subtraction; see details in [START_REF] Girard | Homing in green turtles (Chelonia mydas): oceanic currents act as a constraint rather than as an information source[END_REF]. The swimming movement bouts between two successive water masses-related GNSS relocations are then reconstructed following the standard procedure. Finally, if the geographical information matters (e.g. for relating the animal's location with some habitat characteristics), the fine-scale swimming movement may be converted (through vector addition) into a fine-scale ground-related movement.

Fig. 6. Reconstructed path of a domestic free-ranging horse (based on unpublished data recorded by Jodie Martin). The horse was GPS-tracked for 24 h with a relocation every 10 s, to obtain a faithful path representation (drawn in red). The GPS relocations were subsampled to one every 30 min. to obtain the kind of path representation that is usually obtained when tracking wild animals (drawn in magenta). Then, the data acquired by an on-board tri-axial high-frequency accelero-magnetometer were used to reconstruct the path between 30min GPS relocations (drawn in yellow) using Eqs (19).

When the device is mounted on the collar of terrestrial animal, the device's posture can dramatically differ from the body posture, and vary during movement because the collar is liable to roll around the neck. Fortunately, terrestrial animals can be expected to always move unbanked (λ = 0). Then, provided the device was mounted on the collar with its surge axis orthogonal to the collar plane, which warrants that this axis remains in (or parallel to) the animal's sagittal plane whatever the collar orientation, the body azimuth corresponds to the azimuth of the device and can therefore be computed, similarly to Eq. (5a), as θ = atan 2 ((a

̅ v 2 +a ̅ w 2 ) m ̅ u -a ̅ u (a ̅ v m ̅ v +a ̅ w m ̅ w ), (a ̅ v m ̅ w -a ̅ w m ̅ v ) ||a ̅ ||). (18)
The x and y coordinates on the horizontal plane are then obtained recurrently simply as: arXiv:2310.05820 20

x i = x i-1 + (cos(Θ I -I Θ ) ṽ xi -sin(Θ I -I Θ ) ṽ yi ) L I /Ṽ I (19a) = x i-1 + ((∆X I Σ i ṽ xi +∆Y I Σ i ṽ yi ) ṽ xi -(∆Y I Σ i ṽ xi -∆X I Σ i ṽ yi ) ṽ yi )/(Σ i 2 ṽ xi +Σ i 2 ṽ yi ) y i = y i-1 + (cos(Θ I -I Θ ) ṽ yi + sin(Θ I -I Θ ) ṽ xi ) L I /Ṽ I (19b) = y i-1 + ((∆X I Σ i ṽ xi +∆Y I Σ i ṽ yi ) ṽ yi + (∆Y I Σ i ṽ xi -∆X I Σ i ṽ yi ) ṽ xi )/(Σ i 2 ṽ xi +Σ i 2 ṽ yi ).
with ṽ xi = ṽ i cos(θ i ) and ṽ yi = ṽ i sin(θ i ), L I = ((∆X I ) 2 +(∆Y I ) 2 ) 0.5 , and Ṽ I = (Σ 2 i ṽ xi +Σ 2 i ṽ yi ) 0.5 . An example is shown in Fig. 6. If the z coordinates are relevant (mountainous environment), they can be inferred from the (x,y) coordinates and a Digital Elevation Model. When estimated at a relatively high frequency (e.g. with ∆t = 1 s) the two speed components, ṽ xi and ṽ yi , tends to be very noisy because the rough locomotion mode of terrestrial animals tends to strongly shake the device. Some smoothing of the path is therefore often necessary.

Conclusions

On-board tri-axial high-frequency accelerometers, magnetometers and gyrometers, alone or coupled, have become cheaper and less energy-consuming. This makes it possible to record data on more individuals belonging to smaller species or for longer durations. The present was written with the aim to help the researchers who wish to address questions that require this type of data to compute the key variables and represent their distributions more easily.

Beyond the various trigonometric expressions that are required to compute heading/posture (and possibly relate their changes in terms of body-related rotations) from raw accelero-magnetic data, a recurrent issue with tri-axial high-frequency devices working continuously comes from the huge amount of raw data to manage. It is however worth noting that the successive raw data values acquired at high-frequency are both noisy and highly serially correlated. This means that amount of information they provide is much lower than the amount of data. The derived variables such as heading and posture can be estimated at the same high frequency f as the one used to acquire the raw data using a ∆t-width sliding window which can dramatically smooths the signal (i.e. remove a lot of noise). However, this procedure also dramatically increases the level of serial correlation mechanically (with a sliding window encompassing k = round(f∆t) values, the range of values considered a given time and at the next time shared k-1values), so that the information provided by the change of heading and posture for the very short time 1/f is likely to be insignificant.

Using a ∆t-width jumping rather than sliding window will smooth the signal as well. The amount of data to store and to process in subsequent analyses is then dramatically reduced (by f∆t, as a single value of heading/posture is kept for each time step ∆t), and the level of serial correlation is reduced as well, so that most information present in the signal is preserved. In this case, the time step ∆t should be not only large enough to flatten body oscillations due to limb movements (Shepard et al. 2008a), but also short enough to preserve the highest possible resolution and warrant that (i) changes of heading/posture can be faithfully represented by orthodromic arcs, (ii) roll remains moderate (|ρ|<π) and (iii) reorientation and rolling speeds remains approximately constant (but they can markedly differ between time steps). In-between posture/heading can no longer be directly estimated, but can arXiv:2310.05820 21 be interpolated if required using Eqs (6-8). Furthermore, Eqs. (11-13) highlight how rotations and changes of posture/heading are interrelated for the time step ∆t considered. For technical reasons, posture must be initially computed at high-frequency using a sliding window to assess DBA (at the same high frequency). However, only a single value of posture (through sub-sampling), and of DBA (through averaging) need to be kept for each time step ∆t. It is also worth noting that the activity level can be reliably estimated using a low-frequency accelerometer [START_REF] Benoit | Accelerating across the landscape: the energetic costs of natal dispersal in a large herbivore[END_REF]. A high-frequency tri-axial accelerometer is therefore required only if one also wishes to reliably assess the body posture.

It may be necessary to keep the whole time series of raw data for performing complex analyses such as spectral analyses to detect and characterize regular body oscillations (e.g. [START_REF] Williams | Identification of animal movement patterns using tri-axial magnetometry[END_REF]). However, there often exist simpler and lighter alternatives. For instance, such oscillations in a flying bird, due to wing movements, can be detected and characterized more simply by keeping only two values for each time step ∆t: the mean d̅ and the coefficient of variation, CV, of the delay between successive times at which the heave component of acceleration goes through a balance point equal to -g. Random fluctuations involves that the delay is exponentially distributed (CV=1). Regular oscillations can therefore be detected by CV<<1, and their period can then be estimated as T = 2d̅ (for d̅ >1/f; Nyquist criterion). Based on this approach and on the value of the elevation angle φ, kestrel falcons could be said to hover (CV<0.5, T<150 ms), perform flapped (φ≈0, CV<0.5, T∈ [150,220] ms) or glided (φ≈0, CV>0.5) fly, or to be perched (φ>>0, CV>0.5; García Silveira, unpublished results). For sea turtles, which swim using a forward-backward movement of front flippers, a similar analysis can be done by considering the surge component of acceleration with respect to a balance point equal to sin(φ)g.

As these examples show, some behaviors can be inferred by comparing values derived from the data acquired by bio-loggers to relevant thresholds (see also Shepard et al. 2008b, R. Wilson et al. 2008[START_REF] Fluhr | Space use and time budget in two populations of griffon vultures in contrasting landscapes[END_REF]. However, when quite different behaviors generate relatively similar signals, this simple approach does not work. Machine learning then can help to infer behavior from 'big data' acquired by on-board loggers [START_REF] Jeantet | Behavioural inference from signal processing using animal-borne multi-sensor loggers: a novel solution to extend the knowledge of sea turtle ecology[END_REF][START_REF] Jeantet | Fully convolutional neural network: A solution to infer animal behaviours from multi-sensor data[END_REF]. It should nevertheless be keep in mind that learning will be faster and/or more easily applicable if the machine is fed with biologically relevant parameters such as, for example, DBA and elevation and bank angles rather than with raw acceleration data. Correcting the parameters that apply to the device to obtain those that apply to the body (Eqs 1) should also avoid that the results obtained were blurred by the discrepancy between the two frames of references, which may vary between individuals. In other terms, an approach based on artificial intelligence will work all the better when supported by human intelligence.

Fig. 1 .

 1 Fig. 1. Representation of various axes and angles. Let us consider a bird flying above a body of water. The XY plane, corresponding to the water surface, is horizontal, with the Earth-bound axes X pointing to East and Y pointing to North, while the Earth-bound axis Z points upwards (ENU convention). The animal-bound axes, referred to as the surge, sway and heave axes, correspond to the tail→head, right→left and belly→ back axes, respectively. The 3D orientation is defined by the azimuth, θ, the elevation angle, φ, and the bank angle λ. The azimuth is the horizontal component of heading. It corresponds to the orientation of the surge axis in the XY plane. The elevation angle corresponds to the orientation of the surge axis, measured from the horizontal plane in the vertical plane at azimuth θ. It acts as both the vertical component of heading and the longitudinal component of posture. The bank angle is the transversal component of posture. It is measured from the horizontal level in the animal's frontal plane (i.e. orthogonally to the surge axis) and corresponds to the orientation of the sway axis in this plane. The overall body inclination, ζ, combines the longitudinal (φ) and transversal (λ) components of posture. It corresponds to the deviation of the heave axis from the Z axis. Yaw, α, pitch, β, and roll, ρ, are rotations done around the three main body axes.

  based on the k = round(f∆t) angular speeds, , with the frequency f.

Fig. 2 .

 2 Fig. 2. Change of heading from (θ i-1 , φ i-1 ) to (θ i , φ i ).

  'surface' angles measured in the animal's current frontal plane, i.e. the plane orthogonal to the unit sphere at their respective vertices. In this example of a counter-clockwise reorientation, the values of the different angles are θ i-1 = φ i-1 = π/6, θ i = π/2, φ i = π/3, ω i ≈ π/3.637, δ i i ≈ π/3.256 and δ f i ≈ π/19.023 (replace π by 180 to obtain values in degrees). The product p = cos(φ i-1 ) cos(δ i i ) = cos(φ i )cos(δ f i ) characterizes the inclination of the reorientation plane, φ max = cos -1 (|p)|) ≈ π/2.978, which also corresponds to the elevation angle of the highest point of the great circle to which the orthodromic arc belongs (reached at azimuth θ i-1 + sgn(tan(δ i i )) cos -1 (tan(φ i-1 )/tan(φ max )) = θ i + sgn(tan(δ f i )) cos -1 (tan(φ i )/tan(φ max )) ≈ π/1.784). For a counter-clockwise reorientation, |δ i i | and |δ f i | cannot be larger than φ max while, for a clockwise reorientation, they cannot be smaller than π-φ max , with φ max ∈ [0, π/2]. A key point is that the change of local arc inclination from δ i i to δ f i compels the bank angle to change passively (i.e. without involving any roll) by δ f i -δ i i . If some roll occurs concomitantly, it simply adds to this passive change. Accordingly, the local arc inclination expressed with respect to the coronal plane, which ranges from

Fig. 3 .

 3 Fig. 3. Decomposition of the change of heading in terms of yaw and pitch. This reorientation can be suitably represented as an orthodromic arc drawn on a unit sphere in a pure body-bound frame of reference: the sphere is centered on the animal location, and the three reference axes correspond to the animal's main body axes. The rotation around the heave axis, i.e. yaw (α), and the rotation around the sway axis, i.e. pitch (β), done sequentially in this order, correspond to the two orthogonal components of the orthodromic arc representing the change of heading. They are related to the size (ω) and initial inclination (δ i ) of this arc, and to the bank angle (λ), by the expression [α i •β i ] T = R 1 (δ

  and y i in latitude, ψ i = Ψ I-1 +y i /R. In any case, the I th GNSS-based step (∆X I , ∆Y I , ∆Z I ), with ∆X I = X I -X I-1 , ∆Y I = Y I -Y I-1 , and ∆Z I = Z I -Z I-1 , can be represented by the ground-truth vector L I [Θ I •Φ I ], with length L I = ((∆X I ) 2 +(∆Y I ) 2 +(∆Z I ) 2 ) 0.5 , azimuth Θ I = atan 2 (∆Y I , ∆X I ) and elevation angle Φ I = sin -1 (∆Z I /L I ).

  and therefore |α i | = |β i | for |η i | = π/4 or |η i | = 3π/4 as expected. Otherwise, with roll speed i

  /2] (latitude-like). For this purpose, the horizontal planar grid is replaced by a spherical grid defined in terms of longitude-like and latitude-like value, and the isotropic bivariate Gaussian kernels are replaced by bivariate von Mises-Fisher kernels. In this way, KDE can apply to heading, (θ, φ), change of heading, (δ i , ω), and posture, either (φ, λ) or (ζ, χ), as azimuth, θ, initial arc inclination, δ i , bank angle, λ, and postural orientation, χ, behave as longitudes, elevation angle, φ, behaves as a latitude, and arc size, ω, and overall inclination, ζ, behave as colatitudes (i.e. π/2 -ω and π/2 -ζ behave as latitudes)


Giulia Matilde Ferrante

Satish C Singh

Pranav Audhkhasi

Colin Thomson

Venkata Vaddineni

Milena Marjanović

Supplementary Material for Seismically-derived porosity of deep-sea sediments over the last 74 Ma in the equatorial Atlantic Ocean: Implications for paleoclimate

 Singh, 2019;[START_REF] Audhkhasi | Discovery of distinct lithosphere-asthenosphere boundary and the Gutenberg discontinuity in the Atlantic Ocean[END_REF], (2) the LITHOS-ILAB refraction survey (Profile 1, [START_REF] Grevemeyer | Ocean bottom seismometer and ocean bottom hydrophone seismic refraction and wide-angle data from profile P02 of Maria S. Merian cruise MSM69 with links to sgy data files[END_REF][START_REF] Vaddineni | Evolution of the crustal and upper mantle seismic structure from 0-27 Ma in the equatorial Atlantic Ocean at 2° 43′S[END_REF], and (3) the ILAB-SPARC reflectionrefraction survey (Profiles 2 and 3, [START_REF] Marjanović | Seismic Crustal Structure and Morphotectonic Features Associated With the Chain Fracture Zone and Their Role in the Evolution of the Equatorial Atlantic Region[END_REF][START_REF] Gregory | Serpentinized peridotite versus thick mafic crust at the Romanche oceanic transform fault[END_REF][START_REF] Growe | Seismic structure of the St. Paul fracture zone and late Cretaceous to mid Eocene oceanic crust in the equatorial Atlantic Ocean near 18°W[END_REF]Wang et al., 2022;Wang & Singh, 2022).

During the TransAtlantic ILAB survey, 2750 km of ultra-deep seismic reflection data were acquired onboard WesternGeco M/V Western Trident in 2015, using a 12 km-long multisensor streamer with a receiver spacing of 3.125 m. The energy source consisted of six subarrays, each containing 8 air guns totalling to a volume of 10170 in 3 , deployed at 15 m depth. The shot spacing varied from 50 to 75 m, depending on the age of the lithosphere. Most of the profile lies on the African Plate (0-75 Ma), crosses the ridge-axis, and samples the South American Plate (0-26 Ma). Portions of these MCS data have been used to study the crustal and lithospheric structures along different part of the profile [START_REF] Mehouachi | Water-rich sublithospheric melt channel in the equatorial Atlantic Ocean[END_REF][START_REF] Audhkhasi | Seismic structure of the upper crust from 0-75 Ma in the Equatorial Atlantic Ocean on the African Plate using ultra long offset seismic data[END_REF][START_REF] Marjanović | Seismic Crustal Structure and Morphotectonic Features Associated With the Chain Fracture Zone and Their Role in the Evolution of the Equatorial Atlantic Region[END_REF][START_REF] Audhkhasi | Discovery of distinct lithosphere-asthenosphere boundary and the Gutenberg discontinuity in the Atlantic Ocean[END_REF].

The collocated OBS data were acquired aboard the German R/V Maria S. Merian during the LITHOS-ILAB experiment in 2017. A total of 71 OBSs were deployed along an 1100 km transect, with OBS spacing of 10-20 km. The data were sampled at 250 Hz for all the instruments. The seismic source consisted of six G-gun clusters (12 guns) configured as two sub-arrays with a total volume of 5440 in 3 . The air-gun array was operated at 210 bar, towed at 7.5 m depth and fired every ~400 m along the profile.

The ILAB-SPARC survey was conducted in 2018 onboard the French N/O Pourquoi Pas? and both MCS and OBS data were acquired simultaneously along two perpendicular profiles. The MCS data were acquired using a 6 km-long streamer, towed at 12 m depth, with a receiver spacing of 6.25 m. The energy source consisted of two airgun sub-arrays, each with 9 airguns, with a total volume of 4990 in 3 , towed at 15 m water depth and fired at 300 m interval. Although parts of these profiles were covered during TransAtlanticILAB ultra-deep seismic survey, we decided to use the ILAB-SPARC MCS data to have uniform and coincident coverage with the refraction data. The refraction survey consisted of 85 multi-component OBSs deployed at a spacing of 14-17 km.

S2. VP-VS ratio uncertainty

We calculated the uncertainty associated with where 𝑇 𝑏 -𝑇 𝑠𝑒𝑎 is the seafloor-to-basement two-way travel time (TWTT) beneath the OBS, 𝐷 is the time delay between the P and P-to-S arrivals, dD is the uncertainty associated with the time delay 𝐷 and 𝑑(𝑇 𝑏 -𝑇 𝑠𝑒𝑎 ) is the uncertainty associated with the seafloor-to-basement TWTT sediment thickness 𝑇 𝑏 -𝑇 𝑠𝑒𝑎 .

For conventional travel-time tomography the picking errors on OBS gathers (dD) are normally 20-30 ms but here, because of the strong reverberatory nature of the P-to-S converted waves, we set it to 60 ms.

The uncertainty associated with picking sediment thickness from the MCS data 𝑑(𝑇 𝑏 -𝑇 𝑠𝑒𝑎 ) was calculated as following:

𝑑(𝑇 𝑏 -𝑇 𝑠𝑒𝑎 ) = √(2 * MCS time sampling) 2 + (OBS relocation uncertainty) 2
where the MCS time sampling is equal to 8 ms and the OBS relocation uncertainty is the rms residual error of the least square fitting of water wave arrival picked in each OBS gather.

We found that the uncertainty associated with 𝑉 𝑃 𝑉 𝑆 is between 0.3 and 4.1.

S3. Porosity uncertainty

We calculated the uncertainty associated with porosity propagating the error through Eq. 3 of the main text [START_REF] Lee | Modified Biot-Gassmann Theory for Calculating Elastic Velocities for Unconsolidated and Consolidated Sediments[END_REF]. It depends on the

𝑉 𝑃 𝑉 𝑆
uncertainties and it can be written as follows:

𝑑(𝜙) = 𝐶 1 𝐶 2 𝑛 ( 𝑉 𝑃 𝑉 𝑆 ) -(1+ 1 𝑛 )

𝑑 ( 𝑉 𝑃 𝑉 𝑆 )

where

𝐶 1 = (3𝐾 + 4𝜇) 1 2𝑛 and 𝐶 2 = (3𝜇𝐺 2 ) -1 2𝑛 .
𝐾 and 𝜇 are the bulk and shear modulus of the matrix material and 𝐺 and 𝑛 are constants taking into account the degree of compaction and the effect of clay on velocities, respectively.

Porosity values with relative uncertainty >10% (9 points in total, all of them at crustal ages < 17 Ma) are considered outliers and excluded from the analysis.

S4. Spearman correlation between porosity and climate proxies

We calculated the formal correlation coefficients of the porosity data with 𝛿 18 𝑂, 𝛿 13 𝐶 and 𝐶𝑂 2 climate proxies in the age range 3-37 Ma. To do that, we used the Spearman's rank correlation coefficient, 𝑟 𝑠 , a nonparametric measure of rank correlation, which can assess both linear and non-linear monotonic relationship and it is weakly affected by outliers.

The results show positive correlation between the porosity and 𝐶𝑂 2 (𝑟 𝑠 = 0.665), a negative correlation between the porosity and 𝛿 13 𝐶 (𝑟 𝑠 = -0.445), a weaker but still significant negative correlation between the porosity and 𝛿 18 𝑂 (𝑟 𝑠 = -0.340).

S5. Limit of the method

For this work we have analysed a huge amount of MCS data (~2000 km) and 69 OBSs, covering oceanic crustal age from 3 to 74 Ma. To have consistent results over such a large spatial extent, we had to simplify the calculation. Furthermore, except for ODP sites 662 and 663, which cored down to the Late Pliocene, we do not have any direct measurements of the sediments. Therefore, our analysis required several approximations and assumptions.

First, the formula we used to calculate average

𝑉 𝑃 𝑉 𝑆
in the sediments is valid only in a 1Dapproximation and with the rays travelling vertically within the sedimentary sequence.

However, layering in the sediment creates internal multiples and the long duration of the converted wave can be explained by the effects of wide-angle ringing in sediment wedges arising from dip on the basement-sediment interface. In other words, lateral variations matter and the location of an OBS relative to the underlying sediment pinch-outs or basin flattening can have a significant effect on the delay, strength and duration of the converted wave.

Secondly, we must underline that, in this work, we only considered average 𝑉 𝑃 𝑉 𝑆 and average porosity in the sediments. Yet, over such a large spatial and temporal scale, the sedimentary column will encompass material characterised by a wide range of seismic velocities, density and porosity. We are aware that this is a severe limitation, but we do not have enough information to constrain the variation of such properties with depth. Instead, although we discussed the role of sediments input from Africa, which could be hemipelagic or turbiditic in origin, we made the basic and realistic assumption of dealing only with pelagic sediment. Redistribution of sediment (of the same lithology) via gravity-transport mechanisms, on the other hand, could be particularly important, especially where the seafloor topography is particularly rugged. The topographic highs may have some material removed by slumps and currents, and this material would accumulate in adjacent local basins. Our OBS spacing is not dense enough to appreciate this phenomenon. As the OBS stations were deployed at equal spacing, indistinctly located on topographic highs and lows, both

𝑉 𝑃 𝑉 𝑆
and porosity are estimated on a regional scale. It would be an enormous task to characterise turbidites and slumps along our profiles, but their scientific impact would be limited.

In addition to the above approximations, we have assumed that in the sediments

𝑉 𝑃 𝑉 𝑆
is a function of porosity only. In modelling the porosity, we also assumed the matrix bulk and shear moduli of clay material and we set the dimensional modelling parameters following the modelling studies for unconsolidated sediments with porosities greater than 40% [START_REF] Lee | Modified Biot-Gassmann Theory for Calculating Elastic Velocities for Unconsolidated and Consolidated Sediments[END_REF].

However, we performed a parametric study on these 4 parameters (matrix bulk and shear moduli, n, G), finding that the sensibility of Eq. 3 to these is quite low.

Supplementary Figures

Figure S1: A

Supplementary Tables

S1: Data

Age

  Figure S1: A) Vertical component of seismic data at OBS 08, B) horizontal component. The distance is from the OBS location. The sediment delay is marked. C) MCS image showing the location of OBS 08 along Profile 1.

Figure S2 :

 S2 Figure S2: A) Vertical component of seismic data at OBS 32, B) horizontal component. The distance is from the OBS location. The sediment delay is marked. C) MCS image showing the location of OBS 32 along Profile 1.

Figure S3 :

 S3 Figure S3: A) Vertical component of seismic data at OBS 39, B) horizontal component. The distance is from the OBS location. The sediment delay is marked. C) MCS image showing the location of OBS 39 along Profile 1.

Figure S4 :

 S4 Figure S4: A) Vertical component of seismic data at OBS 59, B) horizontal component. The distance is from the OBS location. The sediment delay is marked. C) MCS image showing the location of OBS 59 along Profile 1.

Figure S5 :

 S5 Figure S5: A) Vertical component of seismic data at OBS 64, B) horizontal component. The distance is from the OBS location. The sediment delay is marked. C) MCS image showing the location of OBS 64 along Profile 2.

Figure S6 :

 S6 Figure S6: A) Vertical component of seismic data at OBS 69, B) horizontal component. The distance is from the OBS location. The sediment delay is marked. C) MCS image showing the location of OBS 69 along Profile 2.

Figure S7 :

 S7 Figure S7: A) Vertical component of seismic data at OBS 75, B) horizontal component. The distance is from the OBS location. The sediment delay is marked. C) MCS image showing the location of OBS 75 along Profile 2.

Figure S8 :

 S8 Figure S8: A) Vertical component of seismic data at OBS 03, B) horizontal component. The distance is from the OBS location. The sediment delay is marked. C) MCS image showing the location of OBS 03 along Profile 3 (modified from Growe et al., 2021).

Figure S9 :

 S9 Figure S9: A) Vertical component of seismic data at OBS 09, B) horizontal component. The distance is from the OBS location. The sediment delay is marked. C) MCS image showing the location of OBS 09 along Profile 3 (modified from Growe et al., 2021).

Figure S10 :

 S10 Figure S10: A) Vertical component of seismic data at OBS 13, B) horizontal component. The distance is from the OBS location. The sediment delay is marked. C) MCS image showing the location of OBS 13 along Profile 3 (modified from Growe et al., 2021).

Contents

Supplementary Text S1-S5 Supplementary Figures S1-S10

Supplementary Tables S1-S2 References Supplementary Text

S1. Seismic surveys and acquisition parameters

The data analysed in this work were collected during three seismic surveys: [START_REF] Grevemeyer | Ocean bottom seismometer and ocean bottom hydrophone seismic refraction and wide-angle data from profile P02 of Maria S. Merian cruise MSM69 with links to sgy data files[END_REF]