Diogo C Barros

Ellen K Longmire

Yi Hui Tee

Motion of finite-size spheres released in a turbulent boundary layer

Keywords: Particle tracking, Particle-laden flow, Turbulent boundary layer

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

Introduction

Particle-laden turbulent flows occur in many applications ranging from industrial processes to the environment such as the pollutant particles in atmosphere, rivers and oceans (e.g., [START_REF] Olsen | Pollutant-particle associations and dynamics in coastal marine environments: a review[END_REF][START_REF] Law | Microplastics in the seas[END_REF][START_REF] Castañeda | Microplastic pollution in St. Lawrence river sediments[END_REF]. In a wallbounded flow, the particle motion is complicated by interactions with both multiple scales in the turbulent flow and the wall itself. The presence of coherent structures such as the alternating high and low momentum regions associated with ejection and sweep events in the near-wall regions [START_REF] Wallace | The wall region in turbulent shear flow[END_REF] can induce complex particle-turbulence interactions. When a particle is larger than the smallest fluid eddies, it can experience variations in shear and normal forces around its circumference. Additionally, wall friction and restitution can affect both the translation and rotation of the particles. Depending on the surrounding conditions, particles can either collide with or lift off from the wall or slide or roll along it. All these effects can significantly impact particle resuspension, deposition and transport which are all very crucial in modeling the particle response in real applications.

Various early experiments examined particle dynamics in turbulent open-channel flows.

Among others, [START_REF] Sutherland | Proposed mechanism for sediment entrainment by turbulent flows[END_REF] investigated how grains in a sediment bed were brought into motion by the fluid. He proposed an entrainment hypothesis whereby strong turbulent eddies could disrupt the viscous sublayer and lift the grain off of the bed. [START_REF] Francis | Experiments on the motion of solitary grains along the bed of a water-stream[END_REF] and [START_REF] Abbott | Saltation and suspension trajectories of solid grains in a water stream[END_REF] observed rolling, saltation and suspension behavior of heavy grains transported over a planar rough bed. Meanwhile, [START_REF] Sumer | Particle motions near the bottom in turbulent flow in an open channel[END_REF] and [START_REF] Sumer | Particle motions near the bottom in turbulent flow in an open channel. Part 2[END_REF] observed that a sand-coated wax sphere with diameter of approximately 30 viscous units (d + ) propagated upwards and downwards repetitively throughout its trajectory in both rough and smooth beds. Hereafter, the symbol + is used to denote quantities normalized by the friction velocity (u τ ) and the kinematic viscosity (ν) of water at 20 • C.

To better understand particle-turbulence interactions in wall-bounded flows, direct visualization techniques were incorporated in investigations by [START_REF] Rashidi | Particle-turbulence interaction in a boundary layer[END_REF], Kaftori et al. c 2020. This manuscript version is made available under the CC-BY-NC-ND 4.0 license. http://creativecommons.org/licenses/by-nc-nd/4.0/ (1995a,b), [START_REF] Niño | Experiments on particle-turbulence interactions in the near-wall region of an open channel flow: implications for sediment transport[END_REF], [START_REF] Van Hout | Spatially and temporally resolved measurements of bead resuspension and saltation in a turbulent water channel flow[END_REF] and [START_REF] Ebrahimian | Dynamics and wall collision of inertial particles in a solid-liquid turbulent channel flow[END_REF] to name a few. These studies concluded that particle resuspension and deposition events in the nearwall regions were strongly influenced by coherent flow structures. Specifically, van Hout (2013) reported that in his time-resolved particle image velocimetry (PIV) and particle tracking velocimetry (PTV) experiments, all lift-off events of polystyrene beads (d + = 10) at friction Reynolds number, Re τ = 435 were due to ejection events generated by passing vortex cores and positive shear. Once lifted beyond the viscous sublayer, the particles either stayed suspended in the fluid or saltated along the wall depending on the type of coherent structures that they encountered. [START_REF] Ebrahimian | Dynamics and wall collision of inertial particles in a solid-liquid turbulent channel flow[END_REF], who studied glass beads (d + = 6.8) in channel flow with Re τ = 410, found the strongest bead accelerations were correlated with fluid ejections and occurred at a wall-normal location of y + = 30. Sweep motions contributed to streamwise decelerations of beads closer to the wall (y + < 20). In addition, [START_REF] Rashidi | Particle-turbulence interaction in a boundary layer[END_REF] and Kaftori et al. (1995a,b), who investigated spanwise motions in dilute particle suspensions in turbulent boundary layers, concluded that the small particles tended to accumulate in the low-speed zones where ejection events were prominent.

A key factor in particle saltation and resuspension is the lift force acting on a particle.

Sphere lift can result from mean shear [START_REF] Saffman | The lift on a small sphere in a slow shear flow[END_REF], individual vortices [START_REF] Auton | The lift force on a spherical body in a rotational flow[END_REF], or turbulence [START_REF] Sutherland | Proposed mechanism for sediment entrainment by turbulent flows[END_REF] in the fluid, from solid body rotation [START_REF] Magnus | Ueber die abweichung der geschosse, und: Ueber eine auffallende erscheinung bei rotirenden körpern[END_REF] or a combination of both (see reviews by [START_REF] Loth | Lift of a solid spherical particle subject to vorticity and/or spin[END_REF][START_REF] Shi | Lift forces on solid spherical particles in unbounded flows[END_REF][START_REF] Shi | Lift forces on solid spherical particles in unbounded flows[END_REF]. When a sphere is close to a wall, the surrounding flow field is distinctly different from that in an unbounded flow. [START_REF] Hall | Measurements of the mean force on a particle near a boundary in turbulent flow[END_REF] and [START_REF] Mollinger | Measurement of the lift force on a particle fixed to the wall in the viscous sublayer of a fully developed turbulent boundary layer[END_REF], who measured the fluid-induced forces acting on a particle fixed to the wall in a turbulent boundary layer using force transducers, reported strong positive mean force in the wall-normal direction.

Hall 's (1988) experimental data showed that, for 3.6 < d + < 140 and particle Reynolds number, 6.5 < Re p < 1250, the normalized mean force could be approximated by F + y = (20.90 ± 1.57)(d + /2) 2.31±0.02 . On the other hand, [START_REF] Zeng | Interactions of a stationary finite-sized particle with wall turbulence[END_REF] fully-resolved direct numerical simulation (DNS) of the flow over a fixed sphere located a distance above the wall showed the opposite. The mean wall-normal forces, evaluated by integrating the nondimensional pressure and viscous stress around the sphere surfaces, were negative for all c 2020. This manuscript version is made available under the CC-BY-NC-ND 4.0 license. http://creativecommons.org/licenses/by-nc-nd/4.0/ spheres with 3.56 ≤ d + ≤ 24.94 centered at a wall-normal location of y + = 17.81 at Re τ = 178.12. Tomographic PIV performed by [START_REF] Van Hout | Experimental study of the flow in the wake of a stationary sphere immersed in a turbulent boundary layer[END_REF] also suggested a negative wall-normal contribution on a tethered sphere with d + = 50 centered at y + = 43 above the wall at Re τ = 352 due to the sphere wake tilting away from the wall. These results imply that the net wall-normal force including the lift force can vary significantly depending on the gap between sphere and wall.

While the presence of a wall can affect the lift force acting on a non-rotating particle considerably, it can also cause the sphere to roll due to frictional torque. In addition, the hydrodynamic torque generated by strong velocity gradients, vorticity and turbulent fluctuations can also cause the sphere to rotate (e.g., [START_REF] Saffman | The lift on a small sphere in a slow shear flow[END_REF][START_REF] Cherukat | A computational study of the inertial lift on a sphere in a linear shear flow field[END_REF][START_REF] Bagchi | Effect of free rotation on the motion of a solid sphere in linear shear flow at moderate Re[END_REF][START_REF] Bluemink | A sphere in a uniformly rotating or shearing flow[END_REF]. In turbulent boundary layer flows, [START_REF] White | Magnus effect in saltation[END_REF] and Niño and García (1994), among others, evaluated the significance of Magnus lift on saltating particles based on comparisons between theoretical and experimental particle trajectories. They concluded that Magnus lift could be a non-negligible part of the overall particle lift force. However, rotation was not quantified in either study.

In many previous numerical simulations of wall-bounded flows, particles were modeled as point-masses with no volume and thus no rotation (e.g., [START_REF] Pedinotti | Direct numerical simulation of particle behaviour in the wall region of turbulent flows in horizontal channels[END_REF][START_REF] Dorgan | Simulation of particles released near the wall in a turbulent boundary layer[END_REF][START_REF] Soldati | Physics and modelling of turbulent particle deposition and entrainment: c 2020[END_REF]. Particle-resolved simulations are relatively limited due to the high computational cost and challenges in getting fine resolution [START_REF] Balachandar | Turbulent dispersed multiphase flow[END_REF]. In most of the particle-laden wall-bounded flow simulations, even though finite-size effects were taken into account, effects due to particle rotation were neither considered nor discussed explicitly (e.g., [START_REF] Pan | Numerical investigation of the effects of large particles on wall-turbulence[END_REF][START_REF] Zeng | Interactions of a stationary finite-sized particle with wall turbulence[END_REF][START_REF] Fornari | The effect of particle density in turbulent channel flow laden with finite size particles in semi-dilute conditions[END_REF]. Several studies including those from [START_REF] Zhao | On particle spin in two-way coupled turbulent channel flow simulations[END_REF], [START_REF] Ardekani | Turbulence modulation in channel flow of finite-size spheroidal particles[END_REF] and [START_REF] Peng | A direct numerical investigation of two-way interactions in a particle-laden turbulent channel flow[END_REF], reported that particle rotation can induce significant effects on the turbulence modulation and should not be neglected.

Several experimental techniques have been proposed to measure sphere translation and rotation simultaneously. For example, [START_REF] Zimmermann | Tracking the dynamics of translation and absolute orientation of a sphere in a turbulent flow[END_REF] and [START_REF] Mathai | Translational and rotational dynamics of a large buoyant sphere in turbulence[END_REF] extracted the sphere position and absolute orientation by comparing the sphere images with unique patterns captured from two perpendicular cameras to a database of synthetic c 2020. This manuscript version is made available under the CC-BY-NC-ND 4.0 license. http://creativecommons.org/licenses/by-nc-nd/4.0/ projections. Meanwhile, [START_REF] Klein | Simultaneous 3D measurement of the translation and rotation of finite-size particles and the flow field in a fully developed turbulent water flow[END_REF] and [START_REF] Barros | Measurement of the translation and rotation of a sphere in fluid flow[END_REF] tracked the markers embedded within or painted over the surface of a sphere using multiple cameras. Then, the angular velocity was computed based on the optimal rotation matrix that best aligned the tracked markers [START_REF] Kabsch | A solution for the best rotation to relate two sets of vectors[END_REF][START_REF] Kabsch | A discussion of the solution for the best rotation to relate two sets of vectors[END_REF].

As presented above, in a particle-laden wall-bounded turbulent flow, the particle-turbulence and particle-wall interactions are complicated to resolve. In the context of particle motion, most attention has been devoted to the two-dimensional translational behavior of particles over short streamwise distances, with little attention to rotation. In order to comprehend the particle dynamics more fully, we track individual spheres over significant streamwise distances while resolving all components of translation and rotation. Multiple sphere densities and flow Reynolds numbers are considered to study the effects of specific gravity and mean shear on sphere motions. Barros et al.'s (2018) methodology is adapted to the requirements of the current experimental setup to reconstruct both the sphere position and orientation. The paper is organized as follows: Section 2 describes the experimental setup, parameters, reconstruction of particle motion and uncertainty analysis; in Section 3, we present the results on sphere translational and rotational kinematics, and discuss how they are affected by both turbulence and the bounding surface; the concluding remarks are summarized in Section 4.

Methodology

Experimental Setup

The experiments were conducted in a recirculating water channel facility at the University of Minnesota. The channel test section, which is constructed of glass, is 8 m long and 1.12 m wide. A 3 mm cylindrical trip-wire was located at the entrance of the test section to trigger the development of a turbulent boundary layer along the bottom wall. Further details on the facility can be found in [START_REF] Gao | Evolution of eddies and packets in turbulent boundary layers[END_REF]. Hereafter x, y and z define the streamwise, wall-normal, and spanwise directions, respectively.

To achieve a repeatable and controllable initial condition, magnetic spheres, with diameter (d) of 6.35 ± 0.05 mm, molded from a mixture of blue machinable wax (913.7 kg m -3 ) and synthetic black iron oxide (5170 kg m -3 ) were used. Here, a small piece of wax was melted inside two hemisphere molds and different amounts of iron oxide particles were added to control the sphere density. Then, the two molds were brought together and chilled immediately to solidify the liquid wax. The resulting spheres were black and opaque. Small markers were painted at arbitrary locations all over the sphere surface using a white oil-based pen for good image contrast (see Fig. 1a). Both the mean inter-marker spacing and mean marker diameter were approximately 0.6 mm.

For each run, a given sphere was held statically on the smooth glass wall in the boundary layer by a magnet at a location 4.2 m downstream of the trip wire and 0.3 m (approximately 4 boundary layer thicknesses, δ) away from the nearest sidewall based on the sphere centroid.

This location will be considered as the origin in x and z, with the bottom wall as y = 0.

A DC 12V 2A push-pull type solenoid from Uxcell was used to hold the cubic neodymium magnet (N40) and the sphere in position. When activated, the solenoid held the plunger such c 2020. This manuscript version is made available under the CC-BY-NC-ND 4.0 license. http://creativecommons.org/licenses/by-nc-nd/4.0/ that the magnet which was connected through a slider will be flush with the outer channel wall. By switching off the power supply, both the plunger and the magnet were retracted with the help of gravity (see Fig. 1b). As the magnet moved away from the bottom wall, the sphere was released, allowing it to propagate with the incoming flow. A coarse screen was located at the end of the test section to capture the sphere and prevent it from recirculating around the channel.

Two pairs of Phantom v210 high-speed cameras from Vision Research Inc. were arranged in stereoscopic configurations to track the sphere in three-dimensional (3D) space over a relatively long field of view. The angle between the two stereoscopic cameras was set to approximately 30 • for both camera pairs. The cameras were positioned with a streamwise overlapping distance of approximately two particle diameters in between their fields of view. All cameras were fitted with 105 mm Nikon Micro-Nikkor lenses with aperture f /16. Scheimpflug mounts were added to all cameras so that the images were uniformly focused across the fields of view. Three white LED panels positioned above the cameras illuminated the domain considered. Prior to running the experiments, the optical system was calibrated by displacing a two-level plate (LaVision Type 22) across nine planes in the spanwise direction for both camera pairs. A third order polynomial fit was obtained for each plane from both cameras pairs using the calibration routine of Davis 8.4 (LaVision GmbH)

to generate the mapping function of the volumetric calibration respectively. The root mean square (r.m.s.) error of the grid point positions was between 0.05 and 0.1 pixels indicating an optimal fit. Image sequences were captured at a sampling frequency of 480 Hz with image resolution of 1280 by 800 pixels, pixel depth of 12-bit, and camera pixel size of 20 µm.

Experimental Parameters

To understand the effect of turbulence, the experiments were conducted at two flow Reynolds numbers. Here, the mean flow statistics of the unperturbed turbulent boundary layers at the initial sphere location were determined from planar PIV measurements in streamwise wall-normal planes. A New Wave Solo II Nd:YAG 532 nm double-pulsed laser system with pulse energy of 30 mJ was used for illumination. The laser sheet illuminated c 2020. This manuscript version is made available under the CC-BY-NC-ND 4.0 license. http://creativecommons.org/licenses/by-nc-nd/4.0/ The PIV images were processed using Davis 7.4 (LaVision GmbH) to obtain the velocity vectors. Normalized cross-correlation [START_REF] Fincham | Low cost, high resolution DPIV for measurement of turbulent fluid flow[END_REF] with an overlap of 50% over initial interrogation window sizes of 64 by 64 pixels followed by three passes of 32 by 32 pixels were employed. All vectors were post-processed with the universal outlier detection criterion [START_REF] Westerweel | Universal outlier detection for PIV data[END_REF] to remove spurious vectors. By applying the Clauser chart method using log-law constant B = 5 and von Kármán constant κ = 0.41, u τ were estimated to be 0.0092 and 0.0193 m s -1 [START_REF] Clauser | The turbulent boundary layer[END_REF][START_REF] Monty | A comparison of turbulent pipe, channel and boundary layer flows[END_REF]. At a water depth of 0.394 m under both fluid conditions, the free-stream velocities (U ∞ ) were Spheres with d = 6.35 ± 0.05 mm and specific gravities (ρ p /ρ f ) of 1.006 (P1), 1.054 (P2) and 1.152 (P3) were considered, where ρ p is the sphere density and ρ f is the fluid density.

This corresponds to d + = du τ /ν of 58 and 122 viscous units, which are comparable to 26 and 46 times the Kolmogorov length scale (η) in the logarithmic region [START_REF] Pope | Turbulent Flows[END_REF]. For each sphere, the density was determined based on the settling velocity (V s ) measured from high-speed imaging of the sphere falling in a quiescent fluid where ρ p = 3C D ρ f V 2 s /4dg + ρ f . Here, C D refers to the drag coefficient obtained from the standard drag curve [START_REF] Clift | Bubbles, drops, and particles[END_REF] while g is the gravitational acceleration. For all spheres, the uncertainty in density was less than 1%. Meanwhile, the initial particle Reynolds numbers defined as Re p = U rel d/ν were 760 and 1840, where U rel = 0.122 and 0.292 m s -1 are the relative velocity between particle and mean fluid at the particle center upon release. In our cases, Stokes numbers (Sk η ), expressed as the ratio of the particle's response time, τ p = (ρ f +2ρ p )d 2 /36νρ f [START_REF] Crowe | Multiphase flow handbook[END_REF] to the characteristic flow time scale based on Kolmogorov, τ η [START_REF] Pope | Turbulent Flows[END_REF], ranged from 59 to 193. Note that although |V s |/U ∞ is relatively small for sphere P1, it is very significant for P2 and P3. Details of the experimental parameters are summarized in Table c 2020. This manuscript version is made available under the CC-BY-NC-ND 4.0 license. http://creativecommons.org/licenses/by-nc-nd/4.0/ 2. The particles were tracked over a streamwise distance up to x ≈ 5.5δ. For each case considered, R = 10 trajectories were captured using the same sphere. For spheres P1 and P3, an additional 55 and 65 runs respectively were completed at both Re τ to track the initial sphere lift-off behavior.

Reconstruction of Particle Translation and Rotation

Before computing the particle translation and rotation, the grayscale images were first pre-processed using Matlab to isolate the sphere from the background. A standard circular Hough Transform routine was applied to locate the sphere. Next, the background surrounding the sphere was removed by setting the intensity values to 0 (black). The extracted sphere images were then imported to Davis 8.4. Here, the images were further processed with 3 x 3 Gaussian smoothing and sharpening filters to increase the dot contrast. Pixel intensity values that were less than those of the white dots were set to 0 to isolate the dots from the sphere image. In all images, the minimum digital dot size was approximately 2 x 2 pixels.

Subsequently, a 3D-PTV routine based on the volumetric calibration mapping function was implemented to reconstruct the dot coordinates from both camera pairs.

The data sets obtained from the PTV were composed of the 3D coordinates of true and c 2020. This manuscript version is made available under the CC-BY-NC-ND 4.0 license. http://creativecommons.org/licenses/by-nc-nd/4.0/ ghost markers and their corresponding 3D velocity vectors. Hence, the filtering methodology proposed by [START_REF] Barros | Measurement of the translation and rotation of a sphere in fluid flow[END_REF] was employed to remove the ghost tracks. Once the true markers had been determined, the sphere centroid was determined by applying the equation of a sphere. Then, a rotation matrix that best aligned the markers of consecutive images was obtained by applying [START_REF] Kabsch | A solution for the best rotation to relate two sets of vectors[END_REF] algorithm. In all runs, at least 8 markers were retained when computing the sphere centroid locations and rotation matrix. Even though all spheres were captured at 480 Hz, different processing frequencies were used in tracking. Depending on the sphere translation speed, the processing frequency was optimized to ensure that marker displacement between images was larger than the disparity uncertainty (see Section 2.4) while limiting maximum sphere streamwise displacement to less than 10 pixels to avoid false marker pairing. More details of the reconstruction process can be found in [START_REF] Barros | Measurement of the translation and rotation of a sphere in fluid flow[END_REF].

Uncertainty Analysis

For the PIV measurements, the random errors of the mean velocity vector in both cases, computed based on a 95% confidence interval, were less than 0.2% of the local mean velocity.

For the r.m.s. of the streamwise velocity fluctuations (u rms ), the corresponding maximum statistical uncertainty estimated based on chi-square analysis with a 95% confidence interval was 3% (see the errorbars plotted in Fig. 2b). The accuracy of δ is limited by the PIV vector spatial resolution and this gives an uncertainty of up to 3% from the estimated value. Meanwhile, the uncertainty of u τ evaluated based on fitting the logarithm law of the wall equation between y + > 30 and y/δ < 0.2 was approximately 1%. These resulted in an uncertainty of 3% in the Re τ calculation. Next, the uncertainty of Re θm , which was estimated based on the statistical error of the mean velocity data, was around 1%.

For the marker reconstruction, the mean disparity error, ǫ * disp calculated by projecting the 3D reconstructed markers back to the camera image in Davis 8.4, was approximately 0.8 px. This gives an estimate of the uncertainty in the marker locations due to reconstruction errors [START_REF] Wieneke | Volume self-calibration for 3D particle image velocimetry[END_REF]. To reduce the noise in computing derivatives, the raw position and orientation data were smoothed by a quintic spline [START_REF] Epps | Evaluating derivatives of experimental data using smoothing splines[END_REF]. The mean uncertain-c 2020. This manuscript version is made available under the CC-BY-NC-ND 4.0 license. http://creativecommons.org/licenses/by-nc-nd/4.0/ ties of the sphere position (x, y and z) as well as orientation (θ x , θ y and θ z ) were computed based on the r.m.s. between the raw and smoothed data [START_REF] Schneiders | Track benchmarking method for uncertainty quantification of particle tracking velocimetry interpolations[END_REF].

For sphere position, the uncertainties were 0.51, 0.18 and 0.95 px; for orientation, the values were 1.2, 1.5 and 0.52 px, respectively. These correspond to mean translational and angular displacement uncertainties of 0.01d and 1.5 • . Lastly, the mean uncertainties of the translational sphere velocities (U p , V p and W p ) were estimated to be 2%, 1% and 4% of the respective U ∞ .

Results and Discussion

Sphere Translation

Sphere Trajectories

Figure 3 shows the sphere wall-normal (y) trajectories plotted against their streamwise distance traveled. In this section, for visual clarity, only results from four out of ten runs (R = 4/10) are plotted for each case, unless otherwise specified. For results from all ten runs, refer to [START_REF] Tee | Three-dimensional tracking of finite-size spheres in a turbulent boundary layer[END_REF]. The upward force required to lift a particle off and away from the wall, which must derive from net pressure and viscous forces, must overcome the constant downward buoyancy force, defined as F b = (ρ p -ρ f )πd 3 g/6, in all cases studied.

In a turbulent boundary layer, the upward force could derive from mean effects such as mean shear or short-term effects related to local fluid motions. Depending on the velocity field surrounding the sphere, wall-normal drag can be significant (see e.g. [START_REF] Van Hout | Spatially and temporally resolved measurements of bead resuspension and saltation in a turbulent water channel flow[END_REF].

Also, if Re p is sufficiently large, i.e. greater than ∼ 100, vortex shedding off of the sphere may generate intermittent wall-normal forces (see [START_REF] Zeng | Interactions of a stationary finite-sized particle with wall turbulence[END_REF][START_REF] Van Hout | Experimental study of the flow in the wake of a stationary sphere immersed in a turbulent boundary layer[END_REF].

Among these possibilities, we first consider the effect of mean shear on the initial sphere wall-normal motion using [START_REF] Hall | Measurements of the mean force on a particle near a boundary in turbulent flow[END_REF] equation obtained for a fixed sphere in a turbulent boundary layer such that F y = F + y ν 2 ρ f . If the mean wall-normal fluid-induced force is greater than the net buoyancy force, then For the least dense sphere P1, at Re τ = 680, F y * = 12. Increasing Re τ to 1320 doubled d + so that the estimated mean force increased sixfold. These estimations correlated very well with our observations for P1 (plotted as black in Fig. 3). At Re τ = 680, this sphere lifted off immediately and within approximately 0.3δ or 3d in 65% and 18% of 75 runs observed, respectively. In the remaining runs, it lifted off at downstream locations significantly greater than 0.3δ. Meanwhile, at Re τ = 1320, this sphere always lifted off immediately upon release in 75 runs observed. Owing to the stronger resultant upward force at higher Re τ , sphere P1 initially rose to greater heights than at the lower Re τ in nine out of ten trajectories reconstructed. Thus, the initial lift-off height correlated strongly with the local mean shear (dU (y)/dy) which is 40% higher for the higher Re τ . In a recent fully-resolved DNS of flow over a sphere initially resting on a porous rough bed, [START_REF] Yousefi | Single sediment dynamics in turbulent flow over a porous bedinsights from interface-resolved simulations[END_REF] showed a resuspension case with F y * = 15 that lifted off initially up to y = 2.5d in one of the four runs simulated.

F y * = (F y -F b )/F b > 0,
Although F y * = -0.2 at Re τ = 680 for sphere P2, it lifted off of the wall one out of ten runs (plotted as green in Fig. 3a). In the remaining runs, this sphere mainly translated along the wall once released. Meanwhile, at higher Re τ (plotted as green in Fig. 3b), this c 2020. This manuscript version is made available under the CC-BY-NC-ND 4.0 license. http://creativecommons.org/licenses/by-nc-nd/4.0/ sphere lifted off of the wall eight out of ten runs due to the stronger mean upward lift force (F y * = 3.4). At this Re τ , all initial lift-off heights were lower than those of P1.

For the densest sphere P3, at the lower Re τ , no initial lift-off was observed in 65 runs.

The sphere did not have sufficient upward lift to overcome the downward buoyance force (F y * = -0.74) and translated along the wall upon release. At the higher Re τ , although F y * = 0.45, sphere P3 lifted off of the wall only once in 65 runs.

In general, [START_REF] Hall | Measurements of the mean force on a particle near a boundary in turbulent flow[END_REF] Spheres that initially lifted off always descended towards the wall after reaching a local maximum in height. While the net buoyancy force was downward in all cases, these descents could be aided by changes in the wall-normal force after the spheres moved away from the wall (see [START_REF] Zeng | Interactions of a stationary finite-sized particle with wall turbulence[END_REF] and additional discussion in Section 3.1.3 below).

Sphere P1 at Re τ = 1320 either recontacted the wall and then lifted off (saltation) or else reached a minimum height above the wall before reascending to a higher location (resuspension). In all other lifting cases, the spheres always recontacted the wall and then either lifted off again (saltation) or else translated along the wall (sliding and/or rolling). These propagation modes were previously observed by [START_REF] Francis | Experiments on the motion of solitary grains along the bed of a water-stream[END_REF], [START_REF] Abbott | Saltation and suspension trajectories of solid grains in a water stream[END_REF], c 2020. This manuscript version is made available under the CC-BY-NC-ND 4.0 license. http://creativecommons.org/licenses/by-nc-nd/4.0/ Sumer and Deigaard (1981), Niño et al. (1994) and van Hout (2013) with particles propagating over rough or smooth beds.

During either saltation or resuspension, sphere P1 frequently ascended to greater heights than those attained after the initial release, with the maximum height observed of y ≈ 2.3d.

This behavior was observed in 10/10 and 7/10 runs investigated at Re τ = 680 and 1320 respectively. At Re τ = 680, this sphere collided with the wall and saltated more frequently than at Re τ = 1320 where it mostly propagated above and further away from the wall. By contrast, the maximum heights observed for spheres P2 and P3 were 0.3d or less throughout their entire trajectories. Interestingly, sphere P3 exhibited a consistent lift-off pattern that was different from the other spheres. Although this sphere did not generally lift off upon release, after it had traveled along the wall a certain distance, repeated small lift-off events began to take place (see figures 5(c, d) for better illustration). This behavior, which was observed for both Re τ , typically began when x 2δ. Even though sphere P1 lifted off to a larger height than other spheres throughout the trajectories, its rising and falling angles were always less than 12 • , and typically stayed below 5 • . For spheres P2 and P3, the rising and falling angles were typically 2 -3 • . and t = 2.262 for 95% confidence level, is always less than ±15% of the respective mean value.

Sphere Streamwise Velocity

Considering spheres that propagated along the wall, the initial acceleration of sphere P3 at Re τ = 680 was much smaller than for sphere P2 at Re τ = 680 and sphere P3 at Re τ = 1320. These trends can be interpreted by taking into account the forward drag and opposing wall friction forces acting on each sphere at the time of release. We estimate the initial drag force using

F D = πρ f C D (U rel d) 2 /8
where C D is the drag coefficient obtained from the standard drag curve [START_REF] Clift | Bubbles, drops, and particles[END_REF] and U rel is the initial relative fluid velocity at the sphere centroid location. Then, the friction force, defined as F f = f c N where f c is a constant friction coefficient and N = F b -F y is the approximate normal force for a non-lifting sphere, is computed for sphere P2 and P3. At Re τ = 680, as the sphere density increased from P2 to P3, friction increased by an order of magnitude. Hence, sphere P3 experienced c 2020. This manuscript version is made available under the CC-BY-NC-ND 4.0 license. http://creativecommons.org/licenses/by-nc-nd/4.0/ stronger opposing wall friction and accelerated more weakly than P2 upon release. As Re τ increased to 1320, the forward drag increased almost fourfold. Due to the smaller forward drag force sphere P3 experienced at Re τ = 680 than at Re τ = 1320, its forward motion was strongly retarded. Note also that for sphere P3, as Re τ increases, the friction force should decrease due to the concomitant increase in F y .

In a uniform, steady, unbounded flow, a sphere would accelerate until it approaches the surrounding fluid velocity. However, in our studies, the presence of both turbulence and the wall modify both the surrounding flow fields and the sphere kinematics. shedding [START_REF] Zeng | Interactions of a stationary finite-sized particle with wall turbulence[END_REF][START_REF] Van Hout | Experimental study of the flow in the wake of a stationary sphere immersed in a turbulent boundary layer[END_REF]. This behavior was observed frequently c 2020. This manuscript version is made available under the CC-BY-NC-ND 4.0 license. http://creativecommons.org/licenses/by-nc-nd/4.0/ for both Re τ , and in almost every run at Re τ = 680 where, downstream of x = 0.5δ, Re p remained larger on average.

Sphere P3 always stayed close to the wall, and where the height varied during very small lift-off events (x > 2δ), it occurs at a much shorter wavelength than the corresponding streamwise velocity fluctuations in the same region (see figures 5(c, d)). Therefore, these variations appear decoupled. When x < 2δ, the velocity curve for Re τ = 680 varied strongly at St ∽ 0.2 also suggesting vortex shedding. This behavior, which was common across all runs for this case, manifests as short wavelengths in Fig. 5c due to the relatively small propagation speed and large relative velocity in that region. After the sphere accelerated again at x ∼ 2δ, the vortex shedding effect was still present. At Re τ = 1320, the streamwise velocity curves also fluctuated at the vortex shedding frequency.

c 2020. This manuscript version is made available under the CC-BY-NC-ND 4.0 license. http://creativecommons.org/licenses/by-nc-nd/4.0/

While the effect of vortex shedding on sphere streamwise velocity is important in most of the cases studied, the simulations by [START_REF] Zeng | Interactions of a stationary finite-sized particle with wall turbulence[END_REF] on a fixed sphere located above the wall showed that the oncoming wall turbulence was also an important contributor to fluctuations in the streamwise force on the sphere. Similarly, in some instances, the timewise variations in velocity magnitude observed herein are likely driven by local high and low momentum regions that surround and move past the sphere. Our preliminary findings on planar fluid velocity fields surrounding a moving sphere at y = 0.75d have shown that a sphere enveloped by a long streamwise region of relatively slow moving fluid traveled significantly slower than the same sphere enveloped by relatively fast moving fluid. Thus, the local fluid velocity variations can explain the strong variations in sphere velocity at a given streamwise location across different runs as well as some occasions where U p /U (y) > 1.

Sphere Wall-normal Velocity

The wall-normal sphere velocities (V p ) at Re τ = 1320 are plotted in Fig. 6a. Two sample V p curves are also plotted in blue on the right axis in Fig. 6b along with their respective wall-normal trajectories in black on the left axis. The results show that upon release, sphere P1 propagated with the highest V p , up to 0.1U ∞ , and hence lifted off to a greater height than other spheres. At the same time, the oscillating frequencies of the wall-normal velocity for P1 in Fig. 6b, as well as in the other runs where Re p were also > 100, matched very well with the Strouhal number associated with vortex shedding. As sphere P2 has greater density than P1, it lifted off with smaller wall-normal velocity. Although the initial V p of P3 was zero, after x ≈ 2δ, the curves began to fluctuate at frequencies corresponding with

St ∼ O(1) or larger, with ascending velocity of approximately 0.05U ∞ . The higher values of St imply that the small repeated lift-off events, which corresponded to the fluctuations in wall-normal velocity, were not prompted by vortex shedding. Surprisingly, even though sphere P3 was the heaviest sphere and lifted off less than 0.25d from the wall, it generally reached V p magnitudes larger than those of P2 and comparable to those of P1. Similar wall-normal velocity trends were observed at the lower Re τ for all spheres (not shown here for brevity). On the other hand, during descents, as shown in Fig. 6a, P1 moved towards the wall with minimum values close to V p ≈ -0.06U ∞ , almost double its settling velocity. In most runs, P1 descended with negative V p value larger than or comparable to those observed for both the denser spheres. A similar finding was observed in the numerical simulation of [START_REF] Yousefi | Single sediment dynamics in turbulent flow over a porous bedinsights from interface-resolved simulations[END_REF], where their resuspending sphere also descended with magnitude larger than V s value at multiple instances along its trajectory.

When the sphere was held at rest in the flow, Re p was 760 and 1840 for Re τ = 680 and 1320, respectively. After strong initial accelerations, these values decreased significantly in all cases. Beyond x = 3δ, the average Re p values were larger for P2 and P3 (O(300) and O(400) at Re τ = 680 and 1320 respectively) than for P1 (O(200) at both Re τ ). While the initial lift-off events were governed mainly by the mean shear lift (as in [START_REF] Hall | Measurements of the mean force on a particle near a boundary in turbulent flow[END_REF] c 2020. This manuscript version is made available under the CC-BY-NC-ND 4.0 license. http://creativecommons.org/licenses/by-nc-nd/4.0/ et al., 2020), the reduction in related upward force associated with the decrease in Re p should weaken the subsequent lift-off activity. However, the resulting wall-normal trajectories and velocities of the spheres demonstrated that the sphere lift-off heights could be larger downstream than those upstream. This signifies the presence of a strong upward impulse of a similar or larger value than that related to the initial mean shear lift. Similarly, for instances where the sphere descended faster than its settling velocity, the presence of a downward fluid impulse can be as significant as that due to the negative buoyancy force.

Yousefi

In this context, [START_REF] Zeng | Interactions of a stationary finite-sized particle with wall turbulence[END_REF] inferred that the strong instantaneous positive and negative wall-normal forces observed on a fixed sphere were closely associated with ejection and sweep events. Studies by [START_REF] Sutherland | Proposed mechanism for sediment entrainment by turbulent flows[END_REF] and van Hout (2013) on moving particles, for example, have shown that ejection events are responsible for particle lift-offs. These findings suggest that the subsequent lift-off events of larger magnitudes observed in the current study are likely triggered and aided by coherent fluid motions. While ejection events could provide the sphere with upward momentum, negative wall-normal fluid motions such as sweeps could provide the sphere with sufficient downward momentum to descend faster than its settling velocity. Depending on the type and strength of coherent structures encountered by a sphere, its wall-normal motion can vary significantly (as shown in Fig. 3) with the change in instantaneous lift and wall-normal drag forces. As a reference, the r.m.s. wall normal fluctuating velocity (v rms ) for 29 < y + < 300 could vary between 3 and 5% of U ∞ at both Re τ .

It is notable that, over the streamwise distances investigated, all sphere trajectories were limited approximately to the buffer and logarithmic regions where turbulent activity is strongest. Thus, although a sphere could be lifted to significant heights above the wall, the wall-normal impulses were not sufficient to drive it to even greater heights. Considering the nature of the impulse, it can derive from a large upward force exerted over a short time or a smaller force exerted over a longer time. For all three spheres examined, the upward impulse must counteract a downward one due to gravity. The wall-normal acceleration for sphere P1 plotted in Fig. 6c revealed multiple behaviors. In one run (solid curves), the sphere acceleration, and net impulse, was initially positive until it reached y + = 150.

c 2020. This manuscript version is made available under the CC-BY-NC-ND 4.0 license. http://creativecommons.org/licenses/by-nc-nd/4.0/ Further downstream, beginning at x = 1.3δ, the sphere again experienced a sustained upward acceleration over the range 175 < y + < 220. By contrast, in the other run (dashed curves), the initial upward acceleration lasted only for a short distance until the sphere rose to y + = 80. Further downstream, however, starting at x = 1.8δ, the same sphere experienced another upward acceleration as it rose from the wall up to y + = 140. The fact that the turbulent activity, and the probability of sustained upward fluid motions, dies off beyond the logarithmic region, helps explain the observed result that the spheres examined never rose above a limit of y + = 282 (y/δ = 0.21) within the field of view investigated.

The wall-normal trajectories plotted in Fig. 3 showed that in most runs, spheres that initially lifted off descended and collided with the wall before lifting off again. Thus, to better understand the effect of collision on sphere lift-off events, the coefficient of restitution (e), which is the ratio of wall-normal velocity after impact (V f ) to wall-normal velocity before impact (V i ), was computed. Among all collision incidents, e was always zero. After colliding with the wall, the spheres slid for a distance of at least 0.1d ( 3 frames) before lifting off again. When the impact Stokes numbers, defined as Sk = ρ p V i d/18νρ f , are calculated, in all trajectories, Sk < 10. Therefore, the e = 0 result for Sk < 10 is consistent with previous experiments on particles moving through quiescent fluid and impacting on solid surfaces [START_REF] Joseph | Particle-wall collisions in a viscous fluid[END_REF][START_REF] Gondret | Bouncing motion of spherical particles in fluids[END_REF].

Sphere Spanwise Motion and Velocity

Sphere trajectories in the streamwise-spanwise plane are plotted in Fig. 7. Once released, instead of propagating along z = 0, the spheres typically moved sideways, occasionally reversing direction, with some crossings over z = 0. For P3 at Re τ = 680 and x < 2δ, the spanwise trajectories exhibited shorter wavelength fluctuations than the other cases. In all cases, the maximum spanwise distance traveled from z = 0 was approximately 12% of the corresponding streamwise distance traveled. This magnitude is similar to that reported by [START_REF] Yousefi | Single sediment dynamics in turbulent flow over a porous bedinsights from interface-resolved simulations[END_REF] for their resuspension case where the sphere traveled mostly away from the bed. Considering all cases, including the lifting sphere P1 at Re τ = 1320, in most runs, the absolute spanwise migration magnitudes were comparable to or larger than the lift-off c 2020. This manuscript version is made available under the CC-BY-NC-ND 4.0 license. http://creativecommons.org/licenses/by-nc-nd/4.0/ trends were observed at the lower Re τ for all spheres (not shown here for brevity). For comparison, at Re τ = 1320 and 29 < y + < 300, the fluid spanwise r.m.s. exceeds the wallnormal r.m.s. value, e.g., at y + = 60, w rms /U ∞ = 0.06 and v rms /U ∞ = 0.04 [START_REF] Jiménez | Turbulent boundary layers and channels at moderate Reynolds numbers[END_REF]. Therefore, within the buffer and logarithmic regions, spanwise forces exerted by the fluid can be larger than the comparable wall-normal forces (see also [START_REF] Zeng | Interactions of a stationary finite-sized particle with wall turbulence[END_REF]. Our preliminary results including fluid velocity fields reveal that the sphere spanwise motion is significantly affected by spanwise fluid drag. For example, we observed that sphere enveloped by the low momentum region tended to travel in the spanwise direction together with the surrounding fluid. More quantitative findings on the correlation between sphere and fluid spanwise velocities are under investigation.

Sphere Rotation

Figure 9 depicts sphere orientation components for R = 4/10 runs of all six cases investigated. The orientations were calculated by integrating the respective angular velocities obtained from the rotation matrix. Depending on the mode of translation, two distinct trends were observed. Spheres that mostly traveled above the wall, namely sphere P1 at both Re τ and sphere P2 at Re τ = 1320, rotated less than half a revolution about all axes throughout the range examined. These weak rotations occurred mostly during the acceleration stage when x < 0.5δ. Notably, even in the presence of strong initial mean shear, the spheres did not develop any significant forward rolling motion (rotation about the negative z-axis). Also, no significant correlations were found between rotation about the streamwise axis (θ x ) and spanwise migration. Without wall friction, these small rotations were most likely induced by the flow structures such as individual vortices or shearing regions.

By contrast, spheres that mostly traveled along the wall, namely sphere P2 at Re τ = 680

and sphere P3 at both Re τ , developed significant rotations starting from x 1.5δ. Upon release, these spheres mainly slid along the wall and barely rolled. on s z = θ x d/2 for a pure rolling motion (plotted as red relative to the left axis) matched very well with the z-trajectories as shown in figures 10 (a, d). This implies that the spanwise sphere motions were induced by the streamwise rotations, and while propagating sideways, these spheres were mostly rolling about the x-axis instead of sliding. These rotations could possibly be initiated by quasi-streamwise vortices or spanwise motions in the fluid.

To quantify the relative importance of sliding and forward rolling, we define a dimensionless rotation rate, α z = |Ω z |d/2U p , which relates the sphere rotational velocity about the spanwise axis, Ω z to its streamwise translational velocity. If α z is 0, the sphere is undergoing pure translation or sliding in the streamwise direction; if α z is 1 and the sphere is in contact with the wall, the sphere is undergoing pure forward rolling without slipping. The dimensionless rotation rates were calculated based on the slope of the mean x-trajectory and θ z curves for the wall-interacting spheres. For the densest sphere P3, when x > 3δ, the mean α z values at Re τ = 680 and 1320 were approximately 0.6 and 0.4 respectively. This c 2020. This manuscript version is made available under the CC-BY-NC-ND 4.0 license. http://creativecommons.org/licenses/by-nc-nd/4.0/ implies a stronger rolling-to-sliding tendency at the downstream location as Re τ decreases.

Although the θ z slope of sphere P2 at Re τ = 680 also increased over time, the mean value of α z ∼ 0.2 infers that translation and sliding remained more important than rotation within the range investigated. Overall, sphere P3 at Re τ = 680 exhibited the strongest rotations.

As discussed in Section 3.1.2, the denser sphere P3 experienced higher wall friction than P2.

Hence, at Re τ = 680, it rolled forward with larger angular velocity than P2. Also, as Re τ increased to 1320, the wall friction weakened and so did the forward rotation.

Finally, we discuss the forward rolling and small repeated lift-off events observed with sphere P3. At Re τ = 680, we noticed that prior to lifting off at x ≈ 2δ, the sphere had

already begun to roll forward and also rotate about both x-and y-axes (see Fig. 10). After the sphere began rolling, the streamwise velocity also increased steeply (see figures 4 and 5c). It is unclear how the rolling initiates. In the sliding P3 case, gravity dominates the wall normal force so that torque from wall friction is relatively constant along the trajectory.

On the other hand, the forward drag force decreases significantly with increasing sphere velocity, and any torque from fluid friction is also likely to change on average. Further, time-wise variations in the approaching fluid velocity field may cause significant variations in torque. Once the sphere begins rolling, Magnus lift can play a role. Kim et al.'s (2014) experiment on spinning spheres in uniform flow reported that as α z increased from 0.1 to 0.6, C L increased from approximately 0.05 to 0.3. These values are comparable to the initial C L estimated from [START_REF] Hall | Measurements of the mean force on a particle near a boundary in turbulent flow[END_REF] equation where C L = 8F y /πρ f (U rel d) 2 ≈ 0.2. Therefore, the contribution of a Magnus force in the P3 lift-off events is likely to be important. It is also notable that the weaker rotations of P2 appeared insufficient to generate these high frequency lift-offs. For P3, the Magnus lift and repeated lift-offs must reduce wall friction and contribute to the second acceleration stage observed.

Conclusions

Three-dimensional translation and rotation of spheres with diameter of 58 and 122 vis- Upon release, the wall-interacting spheres first slid along the wall with minimal rotation.

At Re τ = 680, the initial acceleration of the dense sphere P3 was significantly retarded by the opposing friction force, in contrast to the other spheres that accelerated steeply over a streamwise distance of δ. In this region, vortex shedding led to strong and regular variations in the streamwise velocity. After the wall-interacting spheres propagated downstream by c 2020. This manuscript version is made available under the CC-BY-NC-ND 4.0 license. http://creativecommons.org/licenses/by-nc-nd/4.0/ ≈ 1.5δ, forward rolling (with slipping) followed by repeated small lift-off events of magnitude less than 0.1d began to occur. Since the collisions were completely inelastic (e = 0), and the wall-normal velocity fluctuated at a frequency higher than that expected for vortex shedding, the repeated lift-off events were not caused either by rebounding or vortex shedding. Instead, Magnus lift is most likely the reason behind the lift-off events. While rolling forward, the spheres also rotated about the x-and y-axes, exhibiting coupled rotation behavior.

In all cases, the spheres moved sideways significantly. Due to direct wall interactions, the spanwise motions of the denser spheres were strongly correlated with rotation about the x-axis. Even though the lifting spheres mostly traveled above the wall, they still migrated up to 12% of the streamwise distance traveled. Their spanwise velocities were mostly equal to or larger than the wall-normal velocities and fluctuated at the vortex shedding frequency.

Hence, spanwise forces were substantial in all cases and could be prompted either by rolling about the x-axis, shed vortices, or spanwise fluid motions.

To summarize, the kinematics of spheres in a boundary layer are complicated. The coherent structures in the boundary layer such as high and low momentum regions, as well as vertical and spanwise fluid motions, can have significant effects on both translation and rotation. Since the finite-size spheres propagated toward their own wakes with velocities smaller than the local mean flow, vortex shedding effects were likely significant in all cases investigated such that some of the sphere velocity variations corresponded with expected shedding frequencies. Wall friction was important in impeding acceleration of the denser spheres as well as in prompting the rolling motions. Recently, we have conducted simultaneous stereoscopic PIV and sphere tracking experiments to resolve fluid motions surrounding a moving sphere. These results will be analyzed and discussed in a future communication to complement the current observations and better understand both the particle-wall and particle-turbulence interactions.

c 2020 .

 2020 Fig. 1: Experimental setup. (a) Top view: two pairs of high-speed cameras (S1 and S2) aligned in stereoscopic configurations for capturing the trajectory and rotation of a marked sphere over a long field of view. Inset: example of sphere captured in grayscale with diameter, d spanning 43 pixels or 6.35 mm. (b) Side view: sphere on smooth glass wall in a turbulent boundary layer. Red boxes illustrate the mechanism used in holding (solenoid ON) and releasing (solenoid OFF) the sphere, respectively.

  Fig. 2: Statistics of the unperturbed turbulent boundary layers. (a) Mean streamwise fluid velocity and (b) r.m.s. of streamwise velocity fluctuations. Lines: DNS profiles from Jiménez et al. (2010). Symbols: PIV data. Blue: Re τ = 680; Red: Re τ = 1320. Error bars in (b) indicate the measurement uncertainties.

Fig. 3 :

 3 Fig. 3: Samples of (R = 4/10 runs) sphere wall-normal trajectories at (a) Re τ = 680 and (b) Re τ = 1320 respectively plotted based on the centroid positions. Left axis: plotted in viscous units; Right axis: normalized by the sphere diameter, d. Specific gravity: P1 = 1.006 (black), P2 = 1.054 (green) and P3 = 1.152 (purple).

  equation provided a good estimate as to whether a sphere lifted off initially. The variations in behavior observed for sphere P2 (and the one exception observed for P3) are likely explained by instantaneous variations in streamwise velocity approaching and surrounding them. The r.m.s. streamwise fluid velocities at the initial sphere centroid positions at Re τ = 680 and 1320 were 20% and 15% of the unperturbed local mean value respectively (see Fig.2b). Thus, according to the[START_REF] Hall | Measurements of the mean force on a particle near a boundary in turbulent flow[END_REF] equation, the instantaneous lift forces could vary up to 40% and 26% from F y respectively. Similarly, the instantaneous wall-normal drag force can vary depending on the local velocity fields. These time-dependent variations in the lift and wall-normal drag forces, can also help explain the large variations in initial lift-off heights observed for the lifting spheres released under the same mean fluid condition. Finally, it is important to point out that[START_REF] Hall | Measurements of the mean force on a particle near a boundary in turbulent flow[END_REF] equation was fitted based on experimental data up to Re p = 1250. Hence, our observations on sphere P3 at Re τ = 1320 with Re p = 1840 suggest that[START_REF] Hall | Measurements of the mean force on a particle near a boundary in turbulent flow[END_REF] equation might over-predict the mean wall-normal fluid-induced force for spheres with Re p > 1250.

Figures 4

 4 Figures4(a, b) represent the sphere mean streamwise velocities (U p ) obtained based on averaging U p at each streamwise location over R = 10 curves for each case. Once released, the spheres gained momentum from the hydrodynamic drag force and accelerated rapidly with the flow. At a streamwise location close to 0.5δ, the acceleration magnitudes decreased sharply and the spheres began to approach an approximate mean terminal velocity. From there onward, the mean velocity curves stay relatively flat for all spheres except for sphere P3 at Re τ = 680 where the sphere accelerated again near x ≈ 1.5δ to a new approximate mean terminal velocity. Interestingly, even though the details of sphere motions between P2 and P3 were comparatively different, at both Re τ , their mean velocity curves tend towards similar values. Note that in all cases, after achieving approximate terminal velocity (x > 2δ), the statistical uncertainty in U p , estimated based on ± tσ √ R where σ is the standard deviation

  Based on figures 4(c, d), where U p is normalized by the mean unperturbed fluid velocity at the height of the sphere centroid (U (y)) obtained from PIV, the sphere velocities mostly lagged behind the local mean fluid velocity over the streamwise distance investigated. Sphere P1 came closer to the mean local fluid velocity than the other particles, with mean velocity of approximately 0.9U (y) at both Re τ . For spheres P2 and P3, the curves began to level off around 0.6U (y) and 0.7U (y) at Re τ = 680 and 1320 respectively. These differences are likely due to the fact that P1 is typically detached from the wall while P2 and P3 are not. Moreover, on some occasions, sphere P1 at both Re τ and sphere P2 at Re τ = 1320 traveled faster than the local fluid mean.In all runs, even after attaining an approximate terminal velocity, the individual sphere velocity curves still fluctuate substantially. To better understand the fluctuations in the velocity curves, sample wall-normal trajectories for sphere P1 (black) and P3 (purple) and their corresponding streamwise velocities (blue) at both Re τ are plotted in Fig.5. For sphere P1, the local y-position correlates positively with the respective value of U p in all four runs shown and in the additional runs corresponding to these cases. As sphere P1 ascended, it gained more momentum from the faster moving fluid away from the wall and accelerated; as it descended, it lost momentum due to the slower moving fluid near the wall and thus began to decelerate. When we estimated the oscillating frequencies (f ) of the streamwise velocity curves at both Re τ , in portions of runs where Re p > 100, the resulting Strouhal number, St = f d/U rel corresponded very well with the range of 0.1 to 0.2 associated with vortex

Fig. 5 :

 5 Fig. 5: Sphere wall-normal trajectories, y + (left axis) and streamwise velocities, U p /U ∞ (right axis; blue). (a, c) Re τ = 680; (b, d) Re τ = 1320. Black: P1; Purple: P3. Solid and dashed lines in (a, b) represent two individual runs.

c 2020 .

 2020 Fig. 6: (a) Samples of (R = 4/10 runs) sphere wall-normal velocity at Re τ = 1320, V p normalized by free stream velocity, U ∞ . Black: P1; Green: P2; Purple: P3. (b, c) Sphere P1 at Re τ = 1320; Left axis: wall-normal trajectories, y + (black); Right axis: sphere wall-normal velocity, V p /U ∞ and acceleration, a y,p respectively (blue). Solid and dashed lines represent two individual runs plotted. Red line in (c) represents a y,p = 0.

Fig. 7 :

 7 Fig. 7: Samples of (R = 4/10 runs) sphere spanwise trajectories, z at (a) Re τ = 680 and (b) Re τ = 1320 respectively.

Fig. 9 :

 9 Fig. 9: Samples of (R = 4/10 runs) sphere orientation tracks (θ x , θ y and θ z ) at Re τ = 1320 (a-c) and Re τ = 680 (d-f ) , respectively.

  cous units at Re τ = 680 and 1320 were reconstructed successfully for multiple cases where the relative importance of gravity and mean shear was varied. Particles released rest c 2020. This manuscript version is made available under the CC-BY-NC-ND 4.0 license. http://creativecommons.org/licenses/by-nc-nd/4.0/ were tracked individually over a streamwise distance up to x ≈ 5.5δ. Among all cases, two distinct types of dynamics were observed based on the sphere wall-normal trajectories and rotation behaviors. They included lifting spheres (sphere P1 at both Re τ and sphere P2 at Re τ = 1320) and wall-interacting spheres (sphere P2 at Re τ = 680 and sphere P3 at both Re τ ). Upon release, the lifting spheres separated from the wall in most runs when the mean shear lift force was larger than the net buoyancy force. Due to velocity variations within the turbulent boundary layer, the initial fluid-induced force and lift-off height varied. Since the spheres always reached a maximum height in the logarithmic region before descending toward the wall, it is clear that any initial upward force, possibly induced by mean shear lift, upward moving fluid, or even vortex shedding, decreased significantly as the sphere suspended. After the period of descent, the spheres either ascended again without returning to the wall or else contacted the wall and then slid before lifting off again to heights reaching y ≈ 2.3d (or y + ≈ 135 and 282 at Re τ = 680 and 1320 respectively) for the least dense sphere. Although these spheres mostly traveled above the wall, they nevertheless lagged behind the mean fluid velocity even after attaining an approximate average terminal velocity, propagating with Re p ∼ O(200). Hence, vortex shedding could typically be present, affecting both the surrounding flow fields and the sphere motions. Since the mean shear lift decreases with Re p , all subsequent lift-off events, especially those that reached heights larger than the initial maximum values, are likely prompted by local fluid upwash or local variations in the velocity field, including those related to vortex shedding when Re p > 100. Sometimes the fluid motions acted in the opposite direction, e.g., pushing sphere P1 downward with wallnormal velocity double its gravitational settling velocity. Throughout their propagation, the lifting spheres translated with very weak or minimal rotations about any axis.

Table 1 :

 1 Summary of turbulent boundary layer properties.

	U ∞ [m s -1 ]	u τ [m s -1 ]	δ [mm]	Re τ	Re θm
	0.215	0.0092 ± 0.0001	73 ± 2	680 ± 20	1870 ± 20
	0.488	0.0193 ± 0.0002	69 ± 2	1320 ± 40	3890 ± 40
	through the bottom glass wall had a thickness of 1 mm. The flow was seeded with silver-
	coated hollow glass spheres from Potters Industries LLC with an average diameter and
	density of 13 µm and 1600 kg m -3 , respectively. At both flow Reynolds numbers, 2000
	image pairs were acquired using a TSI Powerview Plus 4MP 16-bit PIV camera with camera
	pixel size of 7.4 µm. The image pairs were captured at a sampling frequency of 1.81 Hz with
	image resolution of 2048 by 2048 pixels.

Table 2 :

 2 Summary of experimental parameters. F y * = (F y -F b )/F b where F y and F b denote the mean wall-normal fluid-induced force based on[START_REF] Hall | Measurements of the mean force on a particle near a boundary in turbulent flow[END_REF] expression and the net buoyancy force, respectively.

	Re τ	U (y = d 2 ) [m s -1 ]	Initial Re p	d +	Sphere	ρ p /ρ f	Sk η |V s |/U ∞	Initial F y	*
	680	0.122	760	58 ± 2	P1	1.006 ± 0.003 59	0.081	12 ± 2
					P2	1.054 ± 0.006 60	0.43	-0.20 ± 0.12
					P3	1.152 ± 0.015 64	0.74	-0.74 ± 0.04
	1320	0.292	1840	122 ± 5	P1	1.006 ± 0.003 176	0.036	72 ± 12
					P2	1.054 ± 0.006 181	0.19	3.4 ± 0.7
					P3	1.152 ± 0.015 193	0.33	0.45 ± 0.24

  and the sphere is expected to lift off upon release from rest on average (see Table2). The uncertainty of this equation was approximately ±0.15F y for both flow Reynolds numbers investigated. c 2020. This manuscript version is made available under the CC-BY-NC-ND 4.0 license. http://creativecommons.org/licenses/by-nc-nd/4.0/

Acknowledgments

The authors thank Nicholas Morse, Ben Hiltbrand and Alessio Gardi for their help with this project. This work was funded by the U.S. National Science Foundation (NSF CBET-c 2020. This manuscript version is made available under the CC-BY-NC-ND 4.0 license. http://creativecommons.org/licenses/by-nc-nd/4.0/ 1510154).


C Fellah

T Douillard

E Maire

S Meille

B Reynard

Gilles Cuny

3D microstructural study of selachimorph enameloid evolution

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

Introduction

Elasmobranchs possess many teeth, frequently replaced, which makes this material quite easy to find and study. Their hard structure allows them to be fossilized more easily than their cartilaginous skeletons, which remain rare in the fossil record [START_REF] Cappetta | Chondrichthyes : Mesozoic and Cenozoic Elasmobranchii : teeth[END_REF].

Enameloid, the hyper-mineralized tissue covering most chondrichthyan teeth, has a complex microstructure. This biomaterial is made from a biological fluorapatite composed of elongated mineral crystallites and a collagen-rich extracellular matrix [START_REF] Legeros | Crystallographic nature of fluoride in enameloids of fish[END_REF][START_REF] Dahm | A Comparative Infrared Spectroscopic Study of Hydroxide and Carbonate Absorption Bands in Spectra of Shark Enameloid, Shark Dentin, and a Geological Apatite[END_REF][START_REF] Daculsi | Ultrastructural study and comparative analysis of fluoride content of enameloid in sea-water and fresh-water sharks[END_REF]. Its mineralization precedes that of the dentine, whereas it is the reverse among tetrapod teeth, which display enamel instead of enameloid. Enamel is produced exclusively by ameloblast cells of epithelial origin, whereas the dentine is the results of odontoblast cells activity, which are of mesenchymal origin. On the other hand, enameloid is attributed to both ameloblasts and odontoblasts [START_REF] Francillon-Vieillot | Microstructure and mineralization of vertebrate skeletal tissues[END_REF][START_REF] Gillis | The homology and phylogeny of chondrichthyan tooth enameloid[END_REF], with a contribution of each of these cells that remains debated. The main component of shark enameloid is fluorapatite (Ca 5 (PO 4 ) 3 F), a fluorine-rich calcium phosphate [START_REF] Baddiel | Spectra structure correlations in hydroxy and fluorapatite[END_REF][START_REF] Bowen | Adhesive Bonding of Various Materials to Hard Tooth Tissues. IV. Bonding to Dentin, Enamel, and Fluorapatite Improved by the Use of a Surface-active Comonomer[END_REF][START_REF] Beevers | The atomic structure of fluor-apatite and its relation to that of tooth and bone material. (With Plates XVI-XVIII[END_REF]. The high fluor content lowers the dissolution potential of this apatite and improves the hardness of the material with respect to hydroxyapatite, with higher stiffness and elastic modulus [START_REF] Gardner | Acoustic microscope study of the elastic properties of fluorapatite and hydroxyapatite, tooth enamel and bone[END_REF]Yoon and Newnham, août, 1969).

Neoselachian teeth exhibit a complex enameloid made of two units (Fig. 1) [START_REF] Enault | Chondrichthyan tooth enameloid: past, present, and future: Chondrichthyan Tooth Enameloid[END_REF][START_REF] Cuny | Evolution of Dental Tissues and Paleobiology in Selachians[END_REF][START_REF] Cuny | The shark fauna from the Middle Triassic (Anisian) of North-Western Nevada[END_REF]. The external unit, named SCE (Single Crystallite Enameloid), is composed of distinct crystallites with a random orientation. Each microcrystal has a 1-2 µm length and a width not exceeding 0.2 µm. The inner unit, named BCE (Bundled Crystallite Enameloid), has crystallites aggregating into bundles, whose diameter is comprised between 5 and 10 µm. Three components are defined according to the organization of the bundles:

-The PBE (Parallel-Bundled Enameloid) where the bundles are oriented parallel to each other and to the surface of the tooth; -The RBE (Radial-Bundled Enameloid) where the bundles are oriented perpendicular to the surface of the tooth; -The TBE (Tangled-Bundled Enameloid) where the bundles intermingle with each other without preferential orientation.

Note that, to reveal ultrastructure, the transverse section was etched 5 s with 4% hydrochloric acid (HCl). The etching produces selective dissolution between the mineralized parts of the tooth and the recrystallized matrix following the fossilization process.

Enameloid microstructure can be studied by scanning electron microscopy (SEM) using two main methods [START_REF] Carter | Effect of etching times and mechanical pretreatment on the enamel of primary teeth: an SEM study[END_REF][START_REF] Seow | The effects of acid-etching on enamel from different clinical variants of amelogenesis imperfecta: an SEM study[END_REF][START_REF] Marshall | SEM Investigation of the Variability of Enamel Surfaces After Simulated Clinical Acid Etching for Pit and Fissure Sealants[END_REF]. The first is to observe the enameloid in longitudinal or transversal section after mechanical cutting of the tooth, followed by polishing and short acid etching of the section. The second is based on successive acid etching in order to map the distribution of the different microstructures over the entire tooth surface. However, these methods permit only a study in two dimensions.

The FIB-SEM tomography (FIB-nt) is widely used in Materials and Life Sciences [START_REF] Giannuzzi | Introduction to Focused Ion Beams: Instrumentation, Theory, Techniques and Practice[END_REF][START_REF] Munroe | The application of focused ion beam microscopy in the material sciences[END_REF]Drobne et al., 2013) but still rarely in Earth Sciences and Palaeontology. This powerful technique has already been used to study the threedimensional structure of dentin [START_REF] Sezen | 3D electron microscopy investigations of human dentin at the micro/nano-scale using focused ion beam based nanostructuring[END_REF][START_REF] Earl | Characterization of dentine structure in three dimensions using FIB-SEM[END_REF][START_REF] Vennat | Three Dimensional Characterization of the Dentin Porous Network Using Confocal Laser Scanning Microscopy[END_REF] and makes it possible to build volumes with a nanometric resolution. The principle of the FIB-nt is based on stacking 2D images, obtained by alternate serial FIB milling and SEM imaging, to reconstruct the microstructure in three dimensions. The volume reconstructed by interpolation between two images can typically reach thousands of µm 3 . Unlike X-ray tomography, inconvenience of FIB-nt is its destructive nature. However, ion milling does not induce selective dissolution as acid etching does.

In this study, the enameloid of an extant carcharhiniform shark was compared to that of an extinct lamniform shark. Although these two sharks are not closely related, previous studies [START_REF] Reif | Morphologie und Ultrastruktur des Hai-"Schmelzes[END_REF] have demonstrated a similar microstructure of their enameloid, allowing to compare and better understand differences due to fossilization processes. These teeth were observed in cross section by conventional scanning electron microscopy in order to obtain initial microstructural information. The FIB-nt technique was applied to obtain 3D reconstructions of the enameloid where SCE meets BCE for both sharks. Consequences for dental microstructure of elasmobranchs hypermineralized tissue are discussed.

Experimental procedure

Materials

Two teeth of modern carcharhiniform sharks, belonging to the genus Carcharhinus, were bought at the Digue market in Antananarivo by one of us (GC) on April 23rd 1996. Their origin was indicated as coming from Toliara, a coastal city situated in the Southwest part of Madagascar. The teeth being isolated, identification at species level was not possible.

The fossil sample, provided by J. Martin, can be referred to Carcharias sp., a lamniform from the Late Cretaceous (Campanian unit), collected at Cap de Naze, Senegal [START_REF] Cuny | A neoselachian shark fauna from the Late Cretaceous of Senegal[END_REF]. The animal tissues used in this article were provided and used with permission and following the ethical guidelines of the University Claude Bernard-Lyon1.

Scanning electron microscopy (SEM)

Classic preparation methods allow observation of the different microstructures of teeth in section or from the surface but is timeconsuming [START_REF] Gillis | The homology and phylogeny of chondrichthyan tooth enameloid[END_REF][START_REF] Enault | Chondrichthyan tooth enameloid: past, present, and future: Chondrichthyan Tooth Enameloid[END_REF]. It is also difficult to evaluate how acid etching alters the global enameloid microstructure. To overcome these problems, the broad ion beam (BIB) milling technique was used to prepare polished surfaces of shark teeth by argon ion bombardment.

The BIB polished sections were prepared as transverse cross-section of the shark teeth apex. Teeth samples were glued with a silver paint onto a titanium blade. The BIB (Gatan Ilion II System) polished cross sections (around 20000 µm 2 ), using an argon beam with an energy set to 6KV for 3 h of milling under vacuum (5.10 -5 Pa).

BIB-polished samples were observed using a FIB-SEM workstation (NVision 40; Carl Zeiss Microscopy) coupling a SIINT zeta ionic column (Seiko Instruments Inc. NanoTechnology, Japan) with a Zeiss Gemini I electronic column. These observations were made in high vacuum, operated at 1 kV with a back-scattered electron (BSE) detector, without coating. SEM images have been processed by ImageJ [START_REF] Collins | ImageJ for microscopy[END_REF], an open-source image processing package, to obtain more microstructural information. The VSNR algorithm [START_REF] Fehrenbach | Variational Algorithms to Remove Stationary Noise: Applications to Microscopy Imaging[END_REF] was used to remove stationary noise and to limit the curtaining effect due to ion abrasion during sample preparation. Subsequently, a thresholding method [START_REF] Xu | Characteristic analysis of Otsu threshold and its applications[END_REF][START_REF] Zhang | Image Segmentation Based on 2D Otsu Method with Histogram Analysis[END_REF] was chosen to create binary images, with elements of interest in white pixels clusters and the background in black, from grayscale images. This method iterates over all possible threshold values and calculates the spread of the pixel levels each side of the threshold to find the optimum.

FIB-SEM tomography

For FIB operation, teeth samples were remounted with a silver paint onto an aluminium stub and coated with a gold layer of about 50 nm by sputtering (EM SCD500, Leica Microsystems GmbH, Germany) to prevent electrostatic charging.

The FIB-SEM serial sectioning experiments were performed with the FIB-SEM workstation (NVision 40; Carl Zeiss Microscopy) used to observe BIB-polished samples. The Atlas 3D software (FIBICS Incorporated, Ottawa, CA), an external scanning system, was used to drive the device.

First, a carbon coating (20 × 20 µm 2 ) was deposited with the insitu gas injection system to define a target area and to protect the surface from ion beam damage. A pattern generated by milling the carbon coating is used for refocusing during acquisition and as markers for image alignment during 3D reconstruction [START_REF] Mainjot | 3D-Characterization of the veneer-zirconia interface using FIB nano-tomography[END_REF]. These marks were enhanced with tungsten deposition and were protected with a last carbon coating. Prior to serial milling and imaging, a coarse trench was milled around the region of interest to a depth of 30 µm.

2D imaging of each slide was conducted using backscattered electrons with the so-called Energy and angle selective BSE detector (EsB) and secondary electrons with the Secondary Electrons Secondary Ions detector (SESI). The FIB-nt experiments were performed at 1.5 kV with a grid potential set to 800 V (EsB). For serial sectioning, FIB milling parameters were adjusted to 30 kV and 1.5nA. These experimental conditions were chosen to fit a cubic voxel size of 20 × 20 × 20 nm 3 . Each stack contains 350 to 700 slices with a target slice thickness of 20 nm.

After raw data acquisition, the image processing was performed in several stages using the ImageJ software. The slices were aligned to correct spatial drifts in the x-and y-directions, due to mechanical and environmental disturbances such as electrostatic charging, with the plugin MultiStackReg. The image alignment results in a loss of voxels at the edges of the volume because of the slice shifting. These useless data were removed to obtain a smooth parallelepiped. The plugin bleach correction was used to remove illumination defect of the stack, such as the shadowing effect. The VSNR algorithm was applied to reduce noise related to curtaining effect on each slice. The treatment involved also other image filters to improve contrast and brightness. After segmentation and thresholding operations, the stack of 2D slices was transformed into a voxel-based data volume where microstructural characteristics were highlighted by labelling. The overall procedure is illustrated in detail in Cougot et al. (Cougot et al., 2018a(Cougot et al., , 2018b)). Other computational operations were directly performed in this volume. Among them, we perfomed Z-projections of the whole stacks of 2D SEM images to visualize 3D data in 2D. We use a method consisting in adding intensity along the z axis for each pixel position. This means superimposing several images in one and having a global view of the details within a volume along an axis.

Results

2D enameloid microstructural organization

In the extant Carcharhinus tooth, contrast is observed between the fluorapatite crystallites and organic matter. The sample preparation method did not reveal the enameloid crystallite clusters but the microstructural pattern of the protein matrix between the crystallites or bundles of the sample (Fig. 2a). In the BCE unit, the crystallite bundles parallel to the tooth surface are separated from each other by protein matrix, represented by the white dots on the segmented image (Fig. 2b). The network of radial bundles surrounded by the protein matrix is also observed. In the upper 2.5 µm of the enameloid (SCE), the bundles of the RBE "open up" like flowers and their crystallites lose their tightly packed, parallel to each other organization. A continuity between crystallites of the RBE and SCE is observed (Fig. 2c). The interest of BIB has been to ensure the conservation of the entire enameloid with respect to etching that may have removed the upper worn and disorganised SCE.

In the fossil Carcharias sample, the three distinct layers usually observed in modern shark teeth (selachimorphs) were not visible as the sample was not acid etched (Fig. 3). The BIB section of the apex did not make it possible to differentiate the different units and microstructural characteristics of the enameloid due to extensive mineralization of the sample after fossilization.

3D microstructural organization of extant Carcharhinus enameloid

In the modern Carcharhinus enameloid, a volume of 18.3 × 17.7 × 6.9 µm 3 volume was acquired in 12 h and a voxel size of 20 × 20 × 20 nm (Fig. 4a). The Z-projection of the 3D representation confirms the results obtained by the 2D microstructural analysis ( § 3.1). The Fig. 4b highlights the continuity of the crystallites network between RBE and SCE units but also the presence of transverse interconnections between radial bundles in the RBE. The microstructural characteristics of Carcharinus enameloid were identified in different colours by labelling operations performed in 3D space (Fig. 4c, 4b and4e). The enameloid structure shows a SCE unit of ~ 1 μm, the upper part is in brown. The spaces between the crystallites oriented parallel to each other and to the surface of the tooth, forming the PBE component of the BCE unit, are shown in cyan. In orange, radial crystallites of the RBE create a periodic lattice and expand to the outermost surface of the tooth. The orthogonal projection just under the surface (Fig. 4f) confirm clearly the network continuity of anastomosed radial bundles in the SCE.

3D microstructural organization of fossil Carcharias enameloid

A 19.0 × 19.1 × 14.0 μm 3 (5080.6 µm 3 ) volume was reconstructed following FIB-SEM tomography experiment on the Carcharias sp. tooth during about 20 h-acquisition time with a similar voxel size (Fig. 5a).

Contrary to the classical SEM images ( § 3.1), SEM images obtained with the FIB-nt routine displays a faint multi-layered structure. In order to enhance this structure, the Z-projection was performed by summing the slices of the 3D volume acquired with SESI and EsB detectors (Fig. 5b). Although the protein matrix has largely decayed in the fossil, the matrix-bundle contrast in the inner part was observed as well as in the upper layer with a 2.5-µm thickness, and displays a similar structure with SCE and BCE as that observed in the extant shark tooth. It was imaged using the same colour code as in Fig. 6 (Fig. 5c, 5d and 5e); RBE is in orange, PBE in cyan and top SCE boundaries in brown. Several cracks in violet (high contrast level on the BSE acquisition) go through the whole sample and are due to geological processes.

The SCE layer (Fig. 5f) is thicker in the fossil Carcharhias teeth than in extant Carcharinus teeth. The difference in SCE thickness between the two samples may have a phylogenetic origin, since the teeth belong to different genera, rather than a diagenetic one. In the volume, radial bundles have many interconnected ramifications with the SCE.

Discussion

FIB-nt allowed obtaining details of the microstructural organization of the enameloid of selachimorphs never observed before. It shows that the structure of the two studied teeth are similar to each other even if the samples belong to two different orders (Carcharhiniformes and Lamniformes) separated by more than 70 million years. This suggests that the structures observed are widely distributed among at least galeomorph sharks and stable through time. The SCE unit appears composed of a wide network of entangled crystallites. This crystallite network is observed to originate from the bundles forming the RBE at the SCE/BCE interface, although it has been suggested that the SCE of selachimorph sharks could be of epithelial origin whereas the BCE was of mesenchymal origin [START_REF] Cuny | The enameloid microstructure of the teeth of Synechodontiform sharks (Chondrichthyes: Neoselachii)[END_REF]. Guinot and Cappetta [START_REF] Guinot | Enameloid microstructure of some Cretaceous Hexanchiformes and Synechodontiformes (Chondrichthyes, Neoselachii): New structures and systematic implications[END_REF] illustrated also a very similar organization in the enameloid of teeth of Hexanchiformes and of Synechodus, although we did not observe the two-layered SCE they describe [START_REF] Guinot | Enameloid microstructure of some Cretaceous Hexanchiformes and Synechodontiformes (Chondrichthyes, Neoselachii): New structures and systematic implications[END_REF]. On the contrary, our results suggest that the whole of the SCE derived from the RBE. The use of FIB methods on Hexanchiformes teeth would certainly allow demonstrating whether the whole of the SCE derives from the RBE in these sharks as well and whether the presence of a double-layered SCE is an artefact due to the limitations of 2D studies. It was recently suggested [START_REF] Delaunois | Microstructural and compositional variation in pacu and piranha teeth related to diet specialization (Teleostei: Serrasalmidae)[END_REF], based on a study of the acrodin of serrasalmid bony fishes, that the cuticle of the acrodin is homologous with the SCE of selachimorph sharks, that the inner enameloid, including the cuticle and dense bundled enameloid, is of epithelial origin while the inner enameloid is of mesenchymal origin. Such a hypothesis would explain the close relationships between SCE and RBE observed here.

Another significant result of the use of SEM-FIB is the clear visualization of frequent branching among the radial bundles forming the RBE, including horizontal link between adjacent bundles. The latter structure was also figured by Guinot and Cappetta (2001) in hexanchiform teeth. Bundles of the RBE are therefore more interconnected than hitherto suspected. This could represent a synapomorphy for the selachimorphs, although more studies will be needed to confirm this hypothesis. Characterizing selachimorphs from their enameloid microstructure has been indeed more problematic since the presence of a PBE in cladodont sharks has been demonstrated by [START_REF] Guinot | Cretaceous stem chondrichthyans survived the end-Permian mass extinction[END_REF].

Conclusion

The 3D reconstruction of enameloid brings a new perspective on this tissue microstructure among selachimorphs. It is first observed that crystallites of the SCE is closely interconnected with those of the RBE although these two components have been proposed to be of different origin. Hence, a more clear-cut boundary between them was hitherto considered. This observation holds for both modern and fossil sharks. Secondly, a 3D reconstruction of the superficial part of the enameloid allows visualizing in detail the branching pattern of the bundles forming the RBE, which is more intricate than hitherto thought. This has the potential to offer a synapomorphy to characterize selachimorph enameloid although more studies will be needed. It should also be investigated whether this branching pattern is also present in the deeper part of the tissue. 3D visualization of enameloid microstructure is therefore a powerful tool for our understanding of the evolutionary pattern of selachimorph enameloid and of the earliest biomineralization processes.

Fig. 1 .

 1 Fig. 1. Tooth of Anomotodon novus from the Ypresian Nanjemoy Formation of Virginia in labial view. Photo by Sten Lennart Jakobsen. The horizontal bar visualizes the transverse section plan (a.). Transverse section of the enameloid of the tooth visualizing the arrangement of the microstructure of the tissue in two units: Single Crystallie Enameloid (SCE), too thin to be properly seen, and Bundled Crystallite Enameloid (BCE). The latter is made of three units: Parallel-Bundled Enameloid (PBE), the bundles of which are oriented perpendicular to the plan of the image, Radial-Bundled Enameloid (RBE) and Tangled-Bundled Enameloid (TBE) (b.). Enlargement of the upper part of the enameloid to better visualize the PBE and RBE (c.).

Fig. 2 .

 2 Fig. 2. Cross-section of a tooth of a modern Carcharhinus after BIB preparation: SEM analysis and image processing. (a.) Overview of the whole cross-section and zoom near the surface in (b.) -(c.) Binary image from thresholding of (b.): Going to the surface, the bundles of the RBE "open up" like flowers and lose their tightly packed, parallel to each other organization. A continuity between crystallites of the RBE and SCE is observed.

Fig. 3 .

 3 Fig. 3. Cross-section of a tooth of a fossil Carcharias after BIB preparation: SEM analysis and image processing. (a.) Overview of the whole cross-section and zoom near the surface in (b.) -(c.) Binary image from thresholding of (b.): Enameloid components of the fossil shark teeth are not highlighted neither by 2D SEM imaging nor by subsequent segmentation.

Fig. 4 .

 4 Fig. 4. 3D FIB/SEM for Modern Carcharhinus enameloid. (a.) Volume after registration of 2D SEM images (18.3 × 17.7 × 6.9 μm 3 , Voxel size: 20 × 20 × 20 nm 3 ) -(b.) Production of a composite 2D image from a Sum -projection in the z direction (i.e. along the axis perpendicular to image plane) of all the slices in the stack. This method allows to clearly make out the enameloid components and highlight transverse few connections between RBE -(c.) Orthogonal view along the Y-axis (after segmentation) of the PBE traces (in cyan) and of the superficial part of the SCE (in brown)-The RBE component is not displayed here. For PBE, the bundles are oriented parallel to each other and to the surface of the tooth -(d.) 3D reconstruction after segmentation of the enameloid components in the volume using contrast thresholding and shape filtering. RBE is in orange, PBE in cyan and top SCE in brown. Note the continuity of the RBE bundles into the SCE area zoomed in (e.) as interconnection between RBE bundles -(f.) Extracted from the volume, left: top view i.e. (along the X-axis) of the segmented black contrasts within SCE (flagged in b.) -right: reconstructed 2D SEM image of the same area in the ZY plane.

Fig. 5 .

 5 Fig. 5. 3D FIB/SEM for fossil Carcharias enameloid (a.) Volume after registration of 2D SEM images (19.0 × 19.1 × 14.0 μm 3 , Voxel size: 20 × 20 × 20 nm 3 ) -(b.) Z-projection to the whole stack of 2D SEM images. The enameloid components are highlighted. The thickness of the SCE is around 2,5µm -(c.) Reconstructed view of both boundary of the SCE and specific contrasts inside (in brown). Several cracks (in violet) go through the whole sample and may be due to geological processes. (d.) 3D reconstruction after segmentation of the enameloid components in the volume using contrast thresholding and shape filtering. RBE is in orange, PBE in cyan and top SCE boundaries in brown. RBE spreads to a great extent into SCE -(e.) Orthogonal view along the Y axis (after segmentation) of the PBE traces (in cyan) and of specific patterns of the SCE (in brown) -The RBE component is not displayed here. PBE keeps contiguous to SCE -(f.) Extracted from the volume, left: top view i.e. (along the X axis) of segmented contrasts within SCE (in brown) -right: reconstructed 2D SEM image of the same area in the ZY plane with entangled contrasts.

Acknowledgements

The authors are very grateful to J. E. Martin (LGL-TPE, CNRS) for providing the tooth sample of Carcharias sp. and the CLyM (Consortium Lyon St-Etienne de Microscopie) for access to the FIB NVision 40. Our two anonymous reviewers' comments allowed us to greatly improve our manuscript.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.