Simon Gsell
email: simon.gsell@univ-amu.fr

Matthias Merkel
email: matthias.merkel@cnrs.fr

Phase separation dynamics in deformable droplets

Phase separation can drive spatial organization of multicomponent mixtures. For instance in developing animal embryos, effective phase separation descriptions have been used to account for the spatial organization of different tissue types. Similarly, separation of different tissue types and the emergence of a polar organization is also observed in cell aggregates mimicking early embryonic axis formation. Here, we describe such aggregates as deformable two-phase fluid droplets, which are suspended in a fluid environment (third phase). Using hybrid finite-volume Lattice-Boltzmann simulations, we numerically explore the out-of-equilibrium routes that can lead to the polar equilibrium state of such a droplet (Janus droplet). We focus on the interplay between spinodal decomposition and advection with hydrodynamic flows driven by interface tensions, which we characterize by a Peclet number P e. Consistent with previous work, for large P e the coarsening process is generally accelerated. However, for intermediate P e we observe long-lived, strongly elongated droplets, where both phases form an alternating stripe pattern. We show that these "croissant" states are close to mechanical equilibrium and coarsen only slowly through diffusive fluxes in an Ostwald-ripening-like process. Finally, we show that a surface tension asymmetry between both droplet phases leads to transient, rotationally symmetric states whose resolution leads to flows reminiscent of Marangoni flows. Our work highlights the importance of advection for the phase separation process in finite, deformable systems.

Introduction

Systems composed of two or more components can undergo phase separation. [1][2][3] During this process, phase domains of different composition first emerge either through spinodal decomposition or nucleation, before these domains coarsen through diffusive processes such as Ostwald ripening. Coarsening is associated with a power-law growth of typical domain size λ, which scales as λ ∼ t 1/3 with time t. 1 In the case of fluid-fluid phase separation, coarsening can also occur through advection with hydrodynamic flows that are driven by the interface tension between adjacent domains. Such hydrodynamic flows can speed up domain coarsening, which can lead to e.g. a linear scaling λ ∼ t. 1 However, the precise interplay between diffusive and hydrodynamic processes can be quite complex, leading to a rich phenomenology. 4 While phase separation has long been employed for industrial processes, it has also been identified as a possible key organizing mechanism in biological systems, from the sub-cellular 5 to the ecological 6 scale. For example, the interior of biological cells is organized in part through fluid compartments that can emerge from phase separation. 5,7,8 Phase separation has moreover been used as effective model to describe the formation of multi-cellular aggregates. 9 But it has also been used to describe separation of different types of biological tissues. [10][11][12][13] This was motivated by experiments showing that in many animal species early differentiated cell populations unmix. 14,15 Phase separation of biological tissues may also help to understand an important question in the field of developmental biology: How does the vertebrate axis form during animal development? To address this and related questions, experimentalists study aggregates of stem cells, called gastruloids, which have been shown to mimic vertebrate development (Fig. 1a,b). [16][17][18] Such gastruloids initially consist of spherical aggregates of undifferentiated stem cells (Fig. 1a). Even under isotropic boundary conditions, these aggregates then reliably 19 form polarized structures with different cell populations (Fig. 1b). [START_REF] Hashmi | bioRxiv[END_REF] This axis later defines the axis of the vertebral column. 17,[START_REF] Van Den Brink | [END_REF] However, how exactly this axis initially emerges has remained unclear.

As a first step towards building a better physical understanding of such cellular aggregates, we study here the passive dynamics of a ternary phase separating system (Fig. 1c-e), where a droplet made out of two fluid phases (1 and 2) is suspended in a third fluid phase (3). Note that such multi-component emulsion droplets [22][23][24] are also important in many other contexts, including the food and pharmaceutical industries. 25 There is a limited number of possibilities for the thermodynamic equilibrium state of such a ternary system, depending on the ratios of the interface tension σ 12 to the surface tensions σ 13 and σ 23 (Fig. 1f). Here, we focus on the case where the equilibrium state is a Janus droplet (Fig. 1d), i.e. where σ 12 < σ 13 + σ 23 and |σ 13 -σ 23 | < σ 12 (white region in Fig. 1f). While the equilibrium state is well known, we are interested in the dynamics that leads to this state, and how it depends on the interplay between diffusion and hydrodynamic flows. We start from an initial state where fluids 1 and 2 are almost homogeneously mixed inside of a spherical droplet (Fig. 1c). Thus, during the phase separation process, the system needs to spontaneously break rotational symmetry to reach the polarized equilibrium state (Fig. 1d).

Here, we numerically study the dynamics of this ternary system using 2D hybrid finite-volume Lattice-Boltzmann simulations of coupled Cahn-Hilliard and Stokes equations. We study the system's dynamics in terms of a Peclet number P e that characterizes the magnitude of hydrodynamic advection as compared to diffusive effects. We find that while increasing P e generally speeds up the phase separation process, intermediate P e give rise to longlived "croissant" states that slow down the relaxation to the final polarized state. These results are not fundamentally changed by including an asymmetry between the surface tensions σ 13 and σ 23 . Taken together, our work demonstrates the importance of hydrodynamic flows in deformable multi-phase droplets.

Model

Free energy

We describe a 2D ternary system where the three components have area fractions C 1 , C 2 , and C 3 , which sum up to one everywhere:

C 1 + C 2 + C 3 = 1. (1)
The phase-separating behavior of these three components is described by the following free energy: 26

F = Ω κ 1 2 C 2 1 (1 -C 1 ) 2 + κ 2 2 C 2 2 (1 -C 2 ) 2 + κ 3 2 C 2 3 (1 -C 3 ) 2 + κ 1 2 (∇C 1 ) 2 + κ 2 2 (∇C 2 ) 2 + κ 3 2 (∇C 3 ) 2 dA.
(

) 2
The integral is over the total system area Ω. For each component j ∈ {1, 2, 3}, the free energy thus includes a quartic double-well potential with local minima at C j = 0 and C j = 1, and an energy density scale κ j . Moreover, gradient terms with coefficients κ j ensure finite widths of all interfaces. For any binary interface between phases of dominating components j and k, i.e. where the third area fraction vanishes, interface width α jk and interface tension σ jk take on the well-known values: 26

α jk = κ j + κ k κ j + κ k and σ jk = α jk 6 (κ j + κ k ). (3)
In the following, we will use the same width α = α jk for all interfaces jk ∈ {12, 13, 23}.

Dynamic equations

To define the dynamics of the area fractions C j such that Eq. ( 1) remains fulfilled, we follow Semprebon et al. 26 by introducing phase fields φ and ψ as

φ = C 1 -C 2 and ψ = C 3 . ( 4
)
Using these definitions together with Eq. ( 1), the free energy in Eq. ( 2) can be expressed in terms of φ and ψ instead of the C j .

We then choose Cahn-Hilliard dynamics for both φ and ψ:

∂φ ∂t + ∇ • (u u uφ) = ∇ • [1 -ψ]Γ φ ∇µ φ (5a) ∂ψ ∂t + ∇ • (u u uψ) = ∇ • Γ ψ ∇µ ψ , ( 5b
)
where t is time, u u u is the local fluid velocity, Γ φ and Γ ψ are mobility parameters, and µ φ = δF/δφ and µ ψ = δF/δψ are the chemical potentials of φ and ψ, respectively. The factor [1 -ψ] in Eq. (5a) ensures that fluids 1 and 2 do not diffuse across fluid 3. This is motivated by the biological background -we exclude any dissociation and re-association of cells (components 1 and 2) from/to the aggregate.

To account for hydrodynamic effects, we describe all phases of the system as incompressible Newtonian fluids, using the incompressible Stokes' equations:

0 = ∇ • u u u, ( 6a
) 0 = -∇Π + ∇ • (2ηS S S) -φ∇µ φ -ψ∇µ ψ . (6b)
Here, Π is the fluid pressure, η is the fluid viscosity, S S S is the shearrate tensor S S S = (∇u u u + (∇u u u) T )/2, and the last two terms correspond to the forces induced by the scalar fields. 1 We focus here on overdamped dynamics, as this is the relevant limit for ∼ 100µmscale cellular aggregates with shear rates on the ∼ day -1 scale.

Because the effective viscosity of embryonic cell aggregates is typically many orders of magnitude above that of water, 27 we choose a ψ-dependent viscosity η(ψ) = (1 -ψ)η i + ψηo, where η i ηo.

Initial and boundary conditions

We initialize the 2D system in a configuration schematized in figure 1c,e. The system initially consists of a circular droplet with diameter L, composed of a mixture of fluids 1 and 2 with ψ = 0, which is surrounded by fluid 3 with ψ = 1 and φ = 0. The inside of the droplet is initialized as φ = φ + ∆φ IC ξ, where φ and ∆φ IC are scalar parameters, and ξ is zero-mean, delta-correlated Gaussian white noise. The parameter φ corresponds to the conserved area average of φ, defined as φ = Ω (1 -ψ)φ dA/ Ω (1 -ψ) dA. We simulate the system in a squared periodic box with side length L t . Parameter values are listed in table 1.

Dimensionless parameters

Equilibrium state

The equilibrium state depends on four dimensionless parameters. These are first, the difference φ in area fractions between fluids 1 and 2. Here we focus exclusively on the case φ = 0, i.e. each fluid occupies half of the droplet area. Second, we set the ratio between droplet size L and interface width α to L/α = 64.

Finally, the equilibrium state also depends on the three tensions σ jk , for which we choose the following two dimensionless ratios:

Ca = 2σ 12 σ 13 + σ 23 , ∆σ = σ 13 -σ 23 σ 12 . ( 7
)
The parameter Ca is akin to a Capillary number -this ratio be- tween the inner and outer surface tensions controls the deformability of the droplet due to capillary effects. Unless stated otherwise, we set Ca = 1. The parameter ∆σ quantifies the surface tension difference between fluids 1 and 2 with respect to fluid 3.

In the following, ∆σ is varied over the range 0 ≤ ∆σ < 1. We are thus in the regime where the equilibrium state is a dipole structure (Fig. 1f).

Relaxation dynamics

The four phenomenological coefficients Γ φ , Γ ψ , η i , and ηo affect only the relaxation dynamics. Absorbing one to rescale time, we obtain three dimensionless parameters. First, we define a Peclet number that compares advective fluxes due to hydrodynamic flows to the diffusive fluxes of the φ-field:

P e = Lα Γ φ η i . ( 8
)
This number is defined on the scale of the droplet size L; it corresponds to the ratio P e = τ D /τ A between a diffusive time scale

τ D = L 2 /D φ with the typical diffusion coefficient scale D φ = Γ φ σ 12 /α ∼ Γ φ (κ 1 + κ 2 )
, and an advective time scale τ A = L/v f with the typical flow velocity scale v f = σ 12 /η i . Increasing the Peclet number corresponds to increasing the influence of advection. The limit case P e = 0 is simulated by ignoring the advection term in equations (5a)-(5b). In this case, we set Γ φ /Γ ψ = 3. In this article, we will vary P e to probe how the interplay of diffusion and hydrodynamic flows affects the phase separation dynamics.

Second, we define a Peclet number with respect to the ψ field, and with respect to the interface width α, as P e ψ = α 2 /(Γ ψ η i ).

Here, we fix this number to a relatively small value of P e ψ = 15/32 to ensure that the outer droplet interface remains stable during the simulations.

Third, we fix the viscosity ratio between droplet and the surrounding fluid to the value η i /ηo = 20. We have checked that changes of this ratio do not strongly affect the simulation results as long as η i /ηo 10.

Numerical implementation

We solve the dynamical equations (5a)-(6b) using a hybrid finitevolume Lattice Boltzmann approach. Similar methods have been used in earlier publications, [28][29][30] using finite-difference schemes for the phase field dynamics. Here, we use a finite-volume scheme instead to ensure exact conservation of the scalar fields φ and ψ. [START_REF] Eymard | Handbook of numerical analysis[END_REF] We discretize space using a Cartesian square grid with spacing ∆x = α. For the advection terms in Eqs. (5a) and (5b), we use a first-order upwind scheme, and for the diffusive fluxes we use a central second-order scheme. We integrate over time using a firstorder explicit Euler scheme. We solve the incompressible Stokes equations (6a) and (6b) by including inertial terms with a small Reynolds number, which we integrate using a two-relaxation-time Lattice Boltzmann method. [START_REF] Gsell | [END_REF] The time step for both finite-volume and Lattice Boltzmann parts is ∆t = (∆x/v f )/640 with the flow velocity scale v f = σ 12 /η i . Details on the numerical method are provided in the Supplementary Information.

Results

In the following, we examine how the droplet polarization dynamics depends (i) on the magnitude of advection as quantified by the Peclet number P e (sections 3.1 and 3.2) and (ii) on surface tension asymmetry, ∆σ = 0 (section 3.3).

Advection speeds up the polarization process

To study the role of advection for the phase separation process, we run simulations with symmetric surface tensions ∆σ = 0 and varying Peclet number P e (Fig. 2, Movies S1-S3). We find well-known features of spinodal decomposition, including an initial decomposition phase where |φ| generally grows until it saturates at ≈ 1, and a subsequent coarsening phase. However, Fig. 2 also suggests (1) that increasing the Peclet number speeds up the phase separation process (discussed below in this section), 1 and (2) that intermediate Peclet numbers give rise to long-lived strongly elongated droplets (discussed in the next section, 3.2).

To quantify in how far increasing the Peclet number speeds up the droplet polarization process, we measure the time evolution of a dipole moment, which we define as Here, c c c is the barycenter of the droplet, defined as c c c = Ω (1ψ)x x x dA/ Ω (1 -ψ) dA. We study the normalized magnitude p = |P P P |/P 0 , where P 0 = L 3 /6 is a reference dipole moment, which corresponds to a circular droplet of diameter L split in two semicircular regions with φ = 1 and φ = -1. We find that indeed, the dipole moment p grows faster for higher Peclet number (Fig. 3a). We moreover define a polarization time tp of the system as the time when the polarization first reaches p ≥ 0.9peq, where peq is the dipole moment of the equilibrium state. Consistent with known results from the literature, 1 we find that advection speeds up the phase separation process (Fig. 3b). However, we note this effect only for P e 100, while tp remains approximately constant for Peclet numbers below 100. To better understand this P e dependence of the polarization process, we now discuss the different phases of the process separately.

P P P (t) = Ω (1 -ψ)φ (x x x -c c c) dA. ( 9
Initial decomposition. We quantify the phase field amplitude |φ| during the initial decomposition phase, which lasts until the time t d ≈ 0.02τ D . To this end, we use the spatial root mean square φ (RMS) of φ, defined as φ2 = Ω (1 -ψ)φ 2 dA/ Ω (1ψ) dA. Fig. 3c shows that φ initially increases exponentially, and reaches a saturation at approximately t d . Except for very large P e, the behavior during this initial phase does not depend very strongly on the Peclet number.

Coarsening. We find that the subsequent coarsening phase is strongly affected by the Peclet number. This is apparent in Fig. 3d, where we plot the time evolution of the typical domain size λ (exact definition in SI). For P e = 0, the domain size λ coarsens according to an approximately constant power law with an exponent close to 1/3. In 2D, this is indeed the predicted exponent for systems that coarsen exclusively through the motion and deformation of domain boundaries, 1 even though we also observe occasional domain evaporation (Fig. 2).

For very large Peclet number, P e 10 3 , we observe a coarsening exponent close to one (Fig. 3d), which corresponds to the expected result for advection-dominated coarsening. For such large P e numbers, advection even induces noticeable coarsening already during the initial decomposition phase (Fig. 3d). This early coarsening is also the likely reason for the slowing down of the initial decomposition for very large P e (Fig. 3c), since diffusive unmixing needs to act over larger length scales in this regime.

The difference in the coarsening behavior between small and large Peclet number is also apparent when studying the behavior of the number of phase domains N D (Fig. 3e). We computationally define a phase domain as a connected region with constant sign of φ. The number of domains N D decreases during the coarsening process until it reaches a minimum of two, which is attained in the equilibrium state. Indeed, N D decreases more rapidly for larger Peclet number Fig. 3e. This confirms that advection speeds up the coarsening process, as N D and λ are strongly correlated (Fig. 3f).

For intermediate Peclet number, we initially observe a relatively fast coarsening (Fig. 3d). However, coarsening subsequently significantly slows down, which is also the reason why the polarization time tp decreases only for very large P e (compare Figs. 3b,d). What could be responsible for this slow down?

Elongated, striped droplets at intermediate P e

We wondered whether the elongated, striped droplets that emerge for intermediate P e (Fig. 2) could be related to the slow down of the coarsening at intermediate P e numbers. To test this hypothesis, we quantified the time evolution of droplet elongation E. We define E by expressing the droplet shape in polar coordinates and computing the second angular Fourier mode (Fig. 4a, details in SI). It is normalized such that E corresponds to the maximal radius variation of the second Fourier mode relative to the initial droplet radius r 0 = L/2.

We find that for P e = 0 and very large P e, the droplet elongation E monotonically increases with time until E ≈ 0.2 (Fig. 4b), which corresponds to the elongation of the final equilibrium state. However, for intermediate P e, droplet elongation first increases as high as E ≈ 0.3 before decreasing back to the final value of ≈ 0.2. Moreover, this regime of transient elongation is correlated in time to the phase of slow coarsening (compare Figs. 3d and4b). This suggests that the elongated droplet structures could indeed be related to the observed coarsening slow down for intermediate P e.

Because the temporal coarsening dynamics differs when varying the P e number, we also plot the droplet elongation over the number of domains N D . This allows to compare droplet elongation E at similar "stages" of the coarsening process across different P e numbers. We find that at the same N D , elongation can be up to two orders of magnitude larger for intermediate P e (Fig. 4c).

We also observe in Fig. 2 that elongated droplets show striped domain patterns. Such a striped arrangement could indeed help to explain the observed coarsening slow down even for intermediate P e. A striped pattern of domains could be close to a mechanical equilibrium where only little hydrodynamic flows occur. In this case, coarsening would occur almost entirely through diffusion in an Ostwald-ripening-like process: Due to an increased Laplace pressure in small stripes as compared to larger stripes of the same color, the larger domain would slowly grow at the expense of the smaller domain through diffusive fluxes across the domains of opposite color.

To test these ideas and to better understand how elongated droplet shapes emerge, we quantify the nematic order of inner domain boundaries (i.e. between droplet domains 1 and 2) by an order parameter Q (Fig. 4d). To this end, we first introduce a tensor field A αβ that characterizes the local orientation of φ interfaces (Fig. 4d):

A αβ = (1 -ψ) ∂ φ ∂xα ∂ φ ∂x β . ( 10
)
Here, Greek letters are dimension indices and we use the nor-

malized field φ = (C 1 -C 2 )/(C 1 + C 2 ) = φ/(1 -ψ).
We then introduce the tensor Q αβ as the spatially averaged, symmetric, traceless part of A αβ :

Q αβ = 1 Aγγ A αβ - 1 2 δ αβ Aγγ . ( 11
)
Here, Einstein notation is used and the brackets • denote spatial averaging over the system. We define the magnitude of the nematic order parameter as

Q = Q 2 xx + Q 2 xy .
When the domain interfaces are randomly oriented, nematic order vanishes, Q = 0, while when all domain interfaces are straight and perfectly aligned, Q = 1.

Generally, the nematic interface alignment Q increases over time, and it does so more quickly for large P e numbers (Fig. 4e). Note that the final equilibrium state has Q ≈ 1, because it consists of a single straight internal interface. Notably, if we plot Q over N D , we find no strong difference across P e numbers (Fig. 4f). Indeed, Fig. 2 also show striped domain patterns not only for intermediate P e, but also for P e = 0. But if the domain pattern is also striped for P e = 0, then why do we find elongated droplets only for intermediate P e (Fig. 4g)?

A possible reason for why P e > 0 droplets are more elongated is that hydrodynamic flows would drive the system closer to mechanical equilibrium, where in mechanical equilibrium the droplets might be elongated. To test these ideas, we first analytically compute the mechanical equilibrium state for a striped droplet. For simplicity, we focus on the limiting case of an infinite chain of equally sized domains (Fig. 4h). In mechanical equilibrium with ∆σ = 0, interfaces between phases 1 and 2 are parallel straight lines, and surfaces to the external fluid 3 are circular arcs. One can furthermore show that the aspect ratio of a single domain in mechanical equilibrium is h/w = 2/Ca (derivation in SI). Thus, in a crude approximation, a droplet with N D domains in mechanical equilibrium would have aspect ratio

N D w/h = N D Ca/2.
To test whether droplets for intermediate P e approach the mechanical equilibrium state, we quantified the domain aspect ratio h/w in our simulations (details in SI). We plot h/w for Ca = 1.5 depending on P e in Fig. 4i and compare to the theoretical solution of the infinite chain (dashed line). Our results suggest that the droplets do indeed approach mechanical equilibrium as P e increases.

To further test these ideas, we compared domain aspect ratios also across different values of Ca (Fig. 4i inset). We found that for smaller Ca, the measured aspect ratio started to become smaller than the theoretical prediction. This makes sense since for smaller Ca, the overall droplet aspect ratio also becomes smaller and so the simulated droplets deviate more strongly from the theoretical picture of a chain of equally sized domains. Indeed, for Ca → 0, the infinite-chain prediction would be h/w = 2/Ca → ∞. However, in this limit, there is no inner interface tension, σ 12 = 0, and the droplet shape is spherical. When dividing a spherical droplet into N D domains of equal width, the domain aspect ratio is h/w ≤ N D , which corresponds to an approximate upper bound for h/w.

Taken together, our results show that for intermediate P e, transient droplet states form, which are striped and elongated. They are close to mechanical equilibrium and coarsen only through diffusive fluxes, which makes them long-lived.

Effect of a difference between the surface tensions of the two droplet phases

We next study the effect of different surface tensions between the two droplet phases, ∆σ = 0. To this end, we carried out simulations with varying 0 ≤ ∆σ < 1 for the Peclet numbers 0, 100, and 10 4 (Fig. 5, Movies S1-S9).

We first examined how varying the surface tension difference ∆σ would change the time tp required to reach the polarized equilibrium state (Fig. 6a). For P e = 0, only simulations with ∆σ = 0 and ∆σ = 0.4 reached the polar equilibrium state within the maximal simulation time of tmax = 50. For P e = 100 and P e = 10 4 , the equilibrium state has always been reached after a finite time tp, no clear trend in the dependency on ∆σ is apparent. To better understand these results, we again separately discuss several phases of the process.

Initial decomposition.

For ∆σ > 0, the initial decomposition into the two droplet phases proceeds from the droplet boundary inwards in concentric circles (Fig. 6b). This phenomenon is called composition wave in the literature: 33,34 For ∆σ > 0, fluid 2 (orange) has a lower surface tension, and so a stripe with an abundance of fluid 2 first appears close to the surface. This supplants fluid 1 (blue) in this region, pushing it further inwards. The slight abundance of fluid 1 next to the outer fluid-2 stripe attracts more of fluid 1, supplanting more of fluid 2, some of which is pushed further inwards, and so on. This initial decomposition proceeds until the whole droplet displays a rotationally symmetric striped pattern, where the stripes reach saturation |φ| ≈ 1 again around

t d ≈ 0.02τ D .
Coarsening. The rotationally symmetric pattern coarsens through progressive broadening of the stripes and occasional annihilation of the innermost phase domain (Figs. 5,6c). For most simulations with P e = 0, the system most often remained in a state with two concentric phase domains until the maximal simulation time (e.g. Fig. 5, for P e = 0, ∆σ = 0.8). However, for finite Peclet number, this rotational symmetry got broken; for P e = 100 during the coarsening phase, and for P e = 10 4 already during the initial decomposition phase. Symmetry breaking. While there is no clear correlation of the polarization time tp with ∆σ, we wondered whether at least the time of symmetry breaking shows a clear dependence on ∆σ. We thus defined the time of symmetry breaking t b as the first dependence on ∆σ. We defined te as the first time when the droplet elongation E reached 90% of its maximum value during any given simulation. This captures both the approximate time of maximum elongation, in cases where a croissant state was created, or the time when E first reached the final plateau value, in cases where the transient flows directly created a 2-phase Janus droplet. We find that te does indeed increase monotonically with ∆σ (Fig. 6e & inset). Thus, while the time tp for the whole process varies nonmonotonically with ∆σ, the time te when the croissant state appears increases monotonically, and is generally much smaller than tp (compare Figs. 6a ande). This means that the variation in tp is due to a variation in the life time of the croissant states.

In the previous section, we concluded that for intermediate P e, these states coarsen only due to an Ostwald-ripening-like process, where smaller domains shrink at the expense of larger ones. This suggests that the variation in tp may just be due to the way the two droplet phases happen to be distributed into different domains by the transient flows that follow the symmetry breaking. Indeed, we find for instance that for P e = 100, the simulations with ∆σ = 0.4 need longest to relax to the 2-domain equilibrium state (Fig. 6a), and these are also the simulations where two approximately equally-sized orange domains form (Fig. 5), leading to an initially slow diffusion between the two domains.

Finally, in simulations with P e = 10 4 , a first symmetry breaking event occurred at a time around 0.01. The ensuing transient flows often resulted in a second approximately rotationally symmetric state with 2 domains at a time around 0.02. This state only later resolved into the polar equilibrium state. We believe this second, 2-domain rotationally symmetric state forms in this case because the first breaking of rotational symmetry occurs before the initial decomposition is finished, i.e. when |φ| is significantly smaller than one within the domains. As a consequence, the interface tension between both domains, which scales as ∼ φ 2 , is smaller, and thus the effective value of ∆σ = (σ 13 -σ 23 )/σ 12 is larger than its nominal value. Hence, we believe that we find the rotationally symmetric 2-domain state as transient state, because the expected equilibrium state is the rotationally symmetric one as long as the effective ∆σ is larger than one (Fig. 1f).

Discussion

In this article, we discussed phase separation in a two-phase deformable droplet. While the thermodynamic equilibrium state for such a system is well known, we studied the route to reach equilibrium. In particular, we focused on the interplay between spinodal decomposition and advection with hydrodynamic flows created by interface tensions. While this interplay has been studied before, 4 it was completely unknown how it plays out in a finite, deformable system. We found that advection can speed up the coarsening process as reported earlier. 1 We moreover found that for intermediate P e numbers, elongated, striped droplet states emerge, which we call "croissant" states. The stripes correspond to a nematic alignment of inner domain boundaries, which emerges for low and intermediate P e. For intermediate P e, the striped droplets elongate as they approach mechanical equilibrium. These states are long-lived, because they coarsen almost exclusively via diffusion, through an Ostwald-ripening-like process. We did not observe such "croissant" states for very large P e numbers, where coarsening is entirely dominated by advection, which leads to coalescence of the phase domains before a striped pattern can develop (Movie S3).

We also studied the effect of an asymmetry between the surface tensions of the two droplet phases, and we found that it changes in particular the beginning of the phase separation process. The emergence of composition waves initially creates a rotationally symmetric system, which then starts to coarsen diffusively. Breaking of this rotationally symmetric state triggers transient flows that are reminiscent of Marangoni flows in droplets. 35,36 These transient flows typically gave rise to a long-lived croissant state.

Taken together, our results demonstrate that advection can play an important role in the coarsening dynamics of deformable twophase droplets.

Symmetry breaking in stem cell aggregates

This work has also important implications for the kinds of biological systems that originally motivated this study. 16,17,[START_REF] Hashmi | bioRxiv[END_REF][START_REF] Van Den Brink | [END_REF][37][38][39] Aggregates of stem cells that mimic early vertebrate axis formation are known to develop a polar structure (Fig. 1a,b), 16 reminiscent of two-phase Janus droplets in thermodynamic equilibrium (Fig. 1d), where the different phases correspond to different cell types. While current studies of such stem cell colonies focus almost exclusively on the biochemical dynamics, 40,41 our work here demonstrates that tissue flow and advection could play an important role as well. Indeed, recent experiments have identified significant large-scale tissue flows in stem cell aggregates. [START_REF] Hashmi | bioRxiv[END_REF] As a first rough estimate, the Peclet number in mouse embryonic stem cell aggregates is on the order of P e ∼ Lσ/Dη ∼ 200. This is based on an aggregate size of L ∼ 200 µm, 20 a tensionto-viscosity ratio of σ/η ∼ 5 µm/h, 42 and a diffusion constant of order D ∼ 5 µm 2 /h (unpublished preliminary data, lab of Pierre-François Lenne). In this P e regime, we observe in our simulations the emergence of the "croissant" states. Indeed, these states appear visually reminiscent of somites, which have also been observed to form in gastruloids, [START_REF] Van Den Brink | [END_REF] while an interface tension has been reported at somite-somite boundaries in vivo. [START_REF] Shelton | bioRxiv[END_REF] However, first, the mechanism of somite formation in vertebrates is generally believed to crucially involve the complex interplay of biochemical signals in a clock-and-wavefront model [START_REF] Oates | [END_REF] (even though purely mechanical models are discussed as well 45 ). Second, different from our model, gastruloids exhibit a polar organization before somites form. [START_REF] Van Den Brink | [END_REF] This raises the question why before polarization there are usually no "croissant" states observed in experiments.

Stem cell cultures are of course much more complex than the model we studied here, and future refined studies should take additional effects into account. First, experiments indicate that the dynamics of some biological signaling molecules and transcription factors is non-conserving in these systems. [START_REF] Hashmi | bioRxiv[END_REF]40,41 To take this into account in our model, a source term needs to be added in the φ equation Eq. (5a). Similarly, we studied here a droplet of constant volume, while the stem cell colonies are typically growing. To implement this, source terms need to be added to the ψ dynamics Eq. (5b) and in the incompressibility condition Eq. (6a).

Second, we focused here on a 2D model, while the stem cell systems motivating this study are 3D aggregates. 16,[START_REF] Hashmi | bioRxiv[END_REF] While there are likely several commonalities between the 2D and 3D cases, we expect differences for instance due to Laplace pressure effects. 34,46,47 Such effects arise only in 3D, where in tube-shaped phase domains, interface tension creates a Laplace pressure that depends on the tube radius. This pressure can lead to necking and breakup of connected domains in a Plateau-Rayleigh instability, 46 but it can also lead to a pumping effect that makes the surface layer of the droplet grow more rapidly. 34 It will be interesting to dissect the consequences of these effects for a deformable 3D droplet in future work.

Third, the cell aggregates are active systems, and so activity may need to be included for an effective description of symmetry breaking in the aggregates. 48,49 Fourth, additional hydrodynamic fields may be required to fully account for some features of the stem cell aggregate dynamics. In particular, the inclusion of more scalar fields in addition to φ may be required, 40,41 or of a polar field as experiments on fish aggregates have pointed to a role of polarity proteins. 50,51 Finally, the real system dynamics is likely influenced by several types of noise. Moreover, a φ-dependence of mechanical properties such as viscosity or local volume growth may be required to explain some of the experimental observations. 52

Fig. 1

 1 Fig. 1 Phase separation and symmetry breaking in deformable multicomponent droplets. (a,b) Microscopy images of cell aggregates mimicking early vertebrate development, taken 1 day apart from each other. Brighter regions indicate the early mesoderm marker T/Brachyury (images taken from Ref. 20 ). (c-e) Ternary fluid mixture model: a droplet, immersed in fluid 3 and initially composed of a mixture of fluids 1 and 2 (c), undergoes phase-separation until reaching the thermodynamic equilibrium state (d) characterized by a surface tension balance between all interfaces. (e) Schematic representation of the physical configuration. (f) Dependence of the thermodynamic equilibrium state on a Capillary number Ca = 2σ 12 /(σ 13 + σ 23 ) and the relative surface tension difference ∆σ = (σ 13 -σ 23 )/σ 12 . The hatched region is unphysical, because at least one of the tensions σ ij would need to be negative there.

Fig. 2

 2 Fig.2Advection affects the coarsening dynamics and the time to reach the polarized equilibrium state. For each value of the Peclet number, a time series of snapshots is shown, representing the φ field from -1 (blue) to 1 (orange). Here we set ∆σ = 0, i.e. the outer surface tension is independent of φ.

Fig. 3

 3 Fig. 3 Advection speeds up droplet polarization by accelerating the coarsening phase. (a) Dipole moment p depending on time for different P e. (b) Polarization time tp as a function of the Peclet number; error bars indicate the standard error of the mean. (c) Spatial root-meansquare of φ depending on time. The black line shows the initial growth with rate ω predicted by linear stability analysis, φ ∼ e ωt . (d) Domain size λ as defined in the SI depending on time. (e) Number of phase domains N D depending on time. (f) Domain size λ as a function of the number of phase domains N D . Legends in panels c-f same as in panel a. In panels c-e, the vertical dashed line indicates the approximate end time t d ≈ 0.02τ D of the initial decomposition phase. The quantities on all ordinate axes are averaged over 25 simulation runs.

4 10 0 10 1 10 2 10 3 10 4 N D = 5 NFig. 4

 4454 Fig. 4 Transient "croissant" states of the droplets, which show an elongated and striped organization of the phase domains. (a) Illustrative example of droplet surface deformation, at P e = 100 and t = 0.1, and resulting angular Fourier spectrum used to define the droplet elongation E. (b) Evolution of the droplet elongation as a function of time. (c) Transient droplet elongation as a function of the Peclet number at fixed numbers of phase domains, N D . (d) Illustrative nematic field based on the φ-field gradient, at P e = 100 and t = 0.1. (e,f) Evolution of the global nematic order as a function of (e) time and (f) N D , for different Peclet numbers. (g) Comparison for two snapshots for P e = 0 and P e = 10 with N D = 5. (h) Limiting case of infinite chain of equally-sized phase domains in mechanical equilibrium, where h/w = 2/Ca (see SI). (i) Phase domain aspect ratio h/w as of function of P e for Ca = 1.5. (i inset) Aspect ratio over Ca for P e = 10. The aspect ratio is averaged for time points with N D = 5 over all domains having exactly two neighboring domains (i.e. all non-branching inner domains). The quantities in all panels are averaged over a series of 25 simulations.

Fig. 6 (

 6 Fig. 6 (a) Polarization time as a function of ∆σ. (b) Composition wave during the initial decomposition phase, for P e = 100 and ∆σ = 0.4. (c) Domain size λ as a function of time. (d) Symmetry breaking time t b over ∆σ. The symmetry breaking time t b is the first time when the nematic order Q of inner domain boundaries increases above a value of 0.01. (d inset) Nematic order Q of the inner domain boundaries over time for different ∆σ. Legend same as in panel c. (e) Time te when a croissant state (or a Janus droplet) emerges over ∆σ. The time te is defined as the first time when droplet elongation E reaches 90% of its maximum value. (e inset) Droplet elongation E over time for different ∆σ. Legend same as in panel c. (f) Directional flow during symmetry breaking. In panels b-f: P e = 100. In each of the panels c-e incl. insets, the respective quantity on the ordinate axis is averaged over 25 simulation runs.

Table 1

 1 Values of dimensionless parameters used in our simulations

	parameter	φ L/α Ca	P e ψ	η i /ηo	∆φ IC	Lt/L
	value	0	64	1	15/32	20	0.01	2
	parameter	∆x/α ∆t/(∆x η i /σ 12 )
	value	1		1/640

[1][2][3][4][5][6][7][8][9][10][11]

1-11 J o u r n a l N a me , [ y e a r ] , [ v o l . ] ,

Acknowledgements

We thank Andrej Košmrlj, Pierre-François Lenne, and Sham Tlili for critical reading of the manuscript. We also thank Pierre-François Lenne and Sham Tlili for useful discussions and experimental data. The project leading to this publication has received funding from the "Investissements d'Avenir" French Government program managed by the French National Research Agency (ANR-16-CONV-0001) and from "Excellence Initiative of Aix-Marseille University -A*MIDEX". The Centre de Calcul Intensif d'Aix-Marseille is acknowledged for granting access to its high performance computing resources.

Surface tension difference

time when the nematic order of the inner domain boundaries, Q, surpassed a value of 0.01. For P e = 100, we found that t b increased monotonically with ∆σ (Fig. 6d & inset).

Transient flows. The symmetry breaking usually triggered complex transient hydrodynamic flows that were qualitatively similar in different parts of the parameter space, and reminiscent of Marangoni flows within a droplet (Fig. 6f). 35,36 For P e = 100, these flows typically resulted in a croissant state.

Slowly relaxing "croissant" state. We wondered whether also the time te it takes for the croissant state to emerge shows a clear

Conflicts of interest

The authors declare no conflict of interest.


Coulouma
email: coulouma@inrae.fr

Assessing new sensor-based volume measurement methods for high-throughput bulk density estimation in the field under various soil conditions

Soil bulk density (BD) is a key soil property in soil science. BD is measured at the scale of the soil horizon with conventional methods : core sampling and rubber balloon. Regardless of the method, BD measurement in the field is cumbersome and time-consuming, especially in stony soils, and new, less invasive methods have emerged but their measurement quality needs to be compared with conventional ones. The photogrammetric technique (SfM) consists of the reconstruction of a given scene in 3D from multi-view photographs. The aim of the present work is to assess the SfM as a rapid, accurate method of measuring the volume of undisturbed soil for BD measurements, regardless of the soil depth and stone content. Ten soil horizons were investigated from various types of soils with different soil properties, especially texture (10 to 48% of clay) and stone content (0 to 73%). The bulk density of each horizon was measured with a reference method (core sampling), an excavation method (rubber balloon) and two recently developed sensor-based methods: SfM and a lightweight flash lidar. For the sensor-based methods, an automated postprocessing method was developed to reduce the human operating time. The BDs measured from the SfM were significantly similar to those measured from core sampling. At this stage, the lightweight flash lidar is not sufficient to measure the BD with high accuracy. Finally, we recommend the use of SfM to measure BD regarding its robustness in varying soil conditions, especially stony soils and we discussed the potential of the method regarding the recent advances on its field use with smartphones.

Introduction

Soil bulk density (BD) is a key soil property in soil science. BD is related to soil porosity, which strongly influences the calculation of carbon content (Goidts et al., 2009), soil available water capacity and nutrient estimation in the soil. Al Shamary et al. (2018)

classify BD

Accepted Article

This article is protected by copyright. All rights reserved.

This article has been accepted for publication and undergone full peer review but has not been through the copyediting, typesetting, pagination and proofreading process which may lead to differences between this version and the Version of Record. Please cite this article as doi: 10.1002/ejss.13115 measurement methods as direct or indirect measurements. The core, excavation and clod methods involve direct measurements of the undisturbed volume of soil samples in the field (core and excavation) or in the laboratory (clod or excavation). Radiation methods such as gamma rays (for example, [START_REF] Bertuzzi | Calibration, field testing, and error analysis of a gamma-ray probe for in situ measurement of dry bulk density[END_REF] rely on indirect measurements in the field and require timeconsuming calibration as well as precautions related to the radiation. Whatever the available method, the measurement of BD is time-consuming and generally requires soil pits to be dug to allow accurate BD measurements. Consequently, BD data are sparse in existing soil databases, which prevent the use of these databases as, for example, inputs for carbon stock calculations [START_REF] Walter | Determining Soil Bulk Density for Carbon Stock calculations: A systematic method comparison[END_REF]. Pedotransfer functions are commonly used for the prediction of BD [START_REF] Qiao | Development of pedotransfer functions for predicting the bulk density in the critical zone on the Loess Plateau, China[END_REF][START_REF] Chen | Building a pedotransfer function for soil bulk density on regional dataset and testing its validity over a larger area[END_REF] based on relationships between well-characterized soil variables (for example, texture and organic carbon content) and bulk density measurements.

The effectiveness of conventional methods used to measure BD at the scale of the soil horizon, such as core sampling and rubber balloon, are also limited by its stone content. The core method is currently considered the standard method for direct BD estimation [START_REF] Blake | Bulk density[END_REF]. However, the reference core volume varies across studies and results in biases in the measured BDs [START_REF] Terry | Effect of soil sample size and included root and wood on bulk density of forested soils[END_REF]. [START_REF] Terry | Effect of soil sample size and included root and wood on bulk density of forested soils[END_REF] found that smaller cores (33 cm 3 ) enhanced the compression effect and that larger cores (2792 cm 3 ) captured macroporosity in the case of forest soils without gravel or stones. Moreover, the core sampling method is not applicable everywhere. The core sampling method is not usable in the case of stony soils or soils without sufficient structure. Excavation methods are then an alternative. The principle of the excavation method is based on the direct measurement in the field of the undisturbed volume with water, sand or Styrofoam balls [START_REF] Grossman | Bulk Density and linear extensibility[END_REF]. To prevent the loss of water due to infiltration, different imaginative ways to seal the border of the excavation have been developed.

The rubber balloon method replaces the void spaces with a rubber balloon filled with a volume of water directly measured on the excavation tool [START_REF] Andraski | Balloon and core sampling for determining Bulk Density of alluvial desert soil[END_REF]). In the same way, Van [START_REF] Van Remortel | Comparison of clod and core methods for determination of Bulk Density[END_REF] sealed the border of the excavation with plastic wrap. Excavation with foam [START_REF] Laundre | Estimating Soil Bulk Density with expanding polyurethane foam[END_REF] or plaster [START_REF] Frisbie | Magnitude and sources of uncertainties in soil organic carbon (SOC) stock assessments at various scales[END_REF][START_REF] Scanlan | Three-dimensional scanning for measurement of bulk density in gravelly soils[END_REF] is also a satisfactory method but is even more time-consuming than the other methods. Moreover, all of the excavation methods recommended in the case of stony soils are difficult to apply vertically [START_REF] Brye | Estimating Bulk Density in vertically exposed stoney alluvium using a modified excavation method[END_REF].

Finally, the clod method consists of sampling clods within the soil horizons and measuring their volume in the laboratory following different methods (Archimedes' principle, X-ray tomography, photogrammetry). The difficulties arise with disaggregation during transport and the risk of selecting specific hard soil structure features within the soil horizon. All of these methods have

Accepted Article

been compared in the literature in terms of their effectiveness and cost as well as the representativeness of the measured BD (for example Al [START_REF] Al-Shammary | Soil Bulk Density estimation methods: A review[END_REF].

New developments have been made in the use of photogrammetry in terms of the number of available platforms, sensors, and software packages [START_REF] Smith | Structure from motion photogrammetry in physical geography[END_REF]. The structure from motion photogrammetric technique, denoted SfM hereafter, consists of the reconstruction of a given scene in 3D from multi-view photographs. [START_REF] Moret-Fernandez | Applicability of the photogrammetry technique to determine the volume and the bulk density of small soil aggregates[END_REF], and more recently [START_REF] Whiting | Rapid bulk density measurement using mobile device photogrammetry[END_REF] characterised the volume of small soil aggregates in the laboratory with SfM techniques . [START_REF] Bauer | A photogrammetric method for calculating soil bulk density[END_REF] tested SfM directly in the field on eight different soil surfaces and compared it with classical sand excavation for stony soils (4 out of 8 sites) and core methods for the rest. The SfM estimated the volume with a better precision than the sand replacement method. The comparison with the core method showed differences, mainly due to soil surface cracks that were not sampled in the cores. In the same way, 3D scanning methods were applied to BD measurement in the laboratory by [START_REF] Rossi | Bulk Density determination by automated three-dimensional laser scanning[END_REF] with a laser scanner on soil aggregates and most recently, [START_REF] Bagnall | An in situ method for quantifying tillage effects on soil structure using multistripe laser triangulation[END_REF] characterised the soil structure from the pit flank with a lidar. [START_REF] Scanlan | Three-dimensional scanning for measurement of bulk density in gravelly soils[END_REF] also used a Kinect sensor from the soil surface but with depth increments to characterize the BD at four different depths within the same excavation. The authors compared the volume calculated from the scanning system with the volume from a classical plaster cast of the excavation hole. The differences were exacerbated in the deeper measurements due to border effects. More recently, [START_REF] Polyakov | An improved excavation method for measuring bulk density of rocky soil using terrestrial LiDAR[END_REF] proposed terrestrial LiDAR as an alternative to the classical excavation method in stony soils for BD measurements at the surface with strong limitations due to the technical requirements and the cost of the device.

However, these promising techniques were only tested in horizontal situations and from the soil surface [START_REF] Bauer | A photogrammetric method for calculating soil bulk density[END_REF][START_REF] Scanlan | Three-dimensional scanning for measurement of bulk density in gravelly soils[END_REF][START_REF] Polyakov | An improved excavation method for measuring bulk density of rocky soil using terrestrial LiDAR[END_REF]. For example, [START_REF] Brye | Estimating Bulk Density in vertically exposed stoney alluvium using a modified excavation method[END_REF] already note the need to measure BD along the exposed face of soil profiles, especially in the case of stony soils. The punctual soil investigations require a vertical approach through the horizonation within the soil pits to obtain a large set of soil property measurements [START_REF] Hartemink | Towards digital soil morphometrics[END_REF]. Currently, there is a lack of study on the comparison of conventional versus new sensor based methods in situ for a large range of soils.

Considering the need for BD data in many soil science applications, the aim of the present work is to assess SfM as an all-terrain method capable of overcoming the limitations of conventional methods for BD measurements, regardless of the stone content. Moreover, a simple and inexpensive 3D sensor is simultaneously tested with photogrammetry and other classical methods used as references: core sampling and the rubber balloon method. To the best of our

Accepted Article

knowledge, our study is the first to assess and compare these four techniques on the same horizons with investigations both from above and from the side of the pit. Starting with a brief description of all the methods tested in the study, their associated protocols and their corresponding operating times, we compared all combinations of methods on the basis of a complete field dataset covering a large range of stone contents, organic carbon, textures, and structures.

Material and methods

Study area

The seven investigated pits that constitute the experimental dataset of this study were located in five study areas, all included in the Langedoc (southern France). They were selected to be representative of the diversity of soil characteristics that may change the bulk density measurement conditions (e.g., coarse fragments, texture, soil depth) and the bulk density itself (e.g., structure). A morphological description was performed on each pit according to the FAO guidelines for soil description [START_REF] Fao | Guidelines for soil description[END_REF]. The final ten investigated horizons were chosen among this available different soil horizons, which is why not all the horizons at the same site were investigated. Five soil pits were dug in January, July and October 2018 for the Pech Rouge (43°08'33'', 3°07'59''), Lavalette (43°08'33'', 3°07'59'') and Mauguio (43°35'33'',4°00'39'') sites, respectively and two additional soil pits were investigated in September 2020 for the Angles (43°31'14'',2°31'45'') and Restinclieres (43°45'58'',3°51'15'') sites. Samples of over 500 g were taken from each horizon for further laboratory characterization. The particle size distribution (texture), gravel and stone content, organic carbon content (OC), and calcium carbonate content (CaCO 3 ) were measured at the French INRAE-ARRAS laboratory. Each sampled horizon is denoted (Site name)_(Pit number when there were multiple pits in a site)_H(horizon number, from top to bottom). Tables 1 and2 summarize the characteristics of this soil dataset. Lavalette_H2, Mauguio_2_H2, and Pech_H2 correspond to deeply tilled horizons with high structural variability.

Lavalette_H4 and Mauguio_1_H3 show medium to high porosity due to high microfaunal activity.

The clay content is highly variable between the investigated horizons, from 10% for Pech_H2 to 48% for Mauguio_3_H2. The organic carbon is generally low except in the Angles_H1 covered by

Accepted Article

This article is protected by copyright. All rights reserved. Douglas trees. Finally, the coarse fragment varied between low content in Mauguio_H1 to a maximum in Mauguio_3_H2. <Table 1 and 2 here>

Soil bulk density estimation

General experimental setup

Each of the ten horizons was investigated both on the side of the pit (core sampling, SfM, flash LiDAR) and from above after digging to reach the given horizon (rubber baloon, SfM, flash LiDAR). <figure 1 here> Investigations from the vertical side of the horizons: core sampling, SfM and Flash LiDAR

The first investigations were performed directly in the side of the pit (figure 1a) with the same procedure repeated for each investigated soil horizon. One unique flash lidar data acquisition of the whole horizon was performed immediately followed by an SfM acquisition that also covered the whole horizon. Both of these procedures were done before soil sampling. Then, five excavations were performed for each sampled horizon. Once the soil was sampled at these five locations, a second flash lidar acquisition of the whole soil horizon and a SfM acquisition of the whole horizon were performed after the soil sampling. Finally, three repetitions of core sampling with 100 cm 3 cores (Blacke and Hartge, 1986) were performed within each of the five excavations.

Hence, the flash lidar data and SfM acquisitions imaged exactly the same volumes, while the core sampling method investigated its own fixed volume samples at the back of each excavation.

Investigations from the horizontal side of the horizons: rubber ballon, SfM and Flash LiDAR

For each investigated horizon, the same procedure described hereafter was used (figure 1b). First, soil was removed to constitute a planar horizontal surface at the average level of the five excavations previously dug from the sides of the pit in order to investigate the same zone (previous section) of a given soil horizon. Five excavations were performed, following the recommendation of the rubber balloon method (NF X31-502). For each excavation, exactly the same excavated volume was measured with three methods: rubber balloon, flash lidar and SfM.

The same procedure was used as above, but with the addition of a density plate around each Accepted Article excavation, fixed with iron fastening clips, which was necessitated by the rubber balloon acquisition (figure 1b). Adding a density plate to the protocol necessitated moving both the density plate and flash lidar from one excavation to another, leading to five acquisitions with each method for a given horizon, compared to only one acquisition from the side of the pit.

The following table (Table 3) summarizes the soil samples taken for each of the ten investigated horizons of the seven different pits. <Table 3 here>

Laboratory sample processing for raw and fine bulk density estimation

The excavated samples were dried (48 h at 105°C) and weighed in the lab. Bulk densities (ρ b ) were determined as the ratio between the dry soil mass and the total sampled volume. As some samples contained high volumes of gravels and stones, we calculated a fine earth bulk density (ρ bFE ) without coarse fragments. To that end, each sample was sieved to extract the coarse fragments (e.g., >2 mm). There were different types of coarse fragments, and their density was measured following the classical method based on the Archimede's principle. The values varied between 1.9 for the highly weathered pebbles in Angles_H1 to 2.6 for the limestones in Restinclieres_H1. The ρ bFE was determined as the ratio between the dry soil mass without the coarse fragments and the total sampling volume without the volume corresponding to the coarse fragments. This latter volume was calculated from the mean bulk density of the coarse fragments given in table 2. Due to the large range of coarse fragments within the dataset, the ρ bFE values were chosen to compare the methods.

Sensor-based volume measurement -data acquisition

For all the investigations, a canopy was used to obtain homogeneous lighting. The canopy was used to prevent shadowing and shade changes during sensor-based data acquisition and soil sampling.

Flash lidar

Accepted Article

We used a PmdTec Pico Flexx flash lidar sensor operated from a smartphone and fixed on an angle iron so that it did not move between the before and after acquisitions. The sensor provides 224 x 171 3D frames with a given precision of 1% of the distance, which corresponds to the order of magnitude of the centimetre in our case (acquisition at approx. one metre).

Acquisition was performed over five seconds for median temporal filtering of the point cloud to decrease noise and eliminate outliers. For a minority of acquisitions, a small offset between the two point clouds was observed, most likely because the angle iron may have shifted during soil sampling. These particular point clouds were thus co-registered by using the ICP algorithm of CloudCompare 2.10.1.

SfM

We used a Nikon D3200 camera with an 18-55 mm AF Nikkor objective set at 18 mm with camera focus fixed so that camera parameters could be considered constant during the acquisitions. A convergent set of at least 20 images -both before and after soil sampling -was acquired to cover the whole area with multi-view stereoscopic imagery, including the scale bars with coded targets.

In the office, all images taken before and after excavation were processed with the Time-SIFT method [START_REF] Feurer | Joining multi-epoch archival aerial images in a single SfM block allows 3-D change detection with almost exclusively image information[END_REF] using Agisoft Photoscan 1.2.6 to produce two 3D point clouds with the exact same geometric reference. Coded targets were used to scale the model. For soil samples taken from the side of the pit, a single pair of before and after point clouds was obtained for the whole horizon. For soil samples taken from above, there was one pair of before and after point clouds for each excavation.

Sensor-based data post-processing

Manual post-processing method

For the manual estimation of excavation volumes, 3D point clouds derived from flash lidar and/or SfM acquisition were processed together in the CloudCompare 2.10.1 software. A cloud-to-cloud distance was computed, and the 3D point clouds were coloured with this distance information.

Each excavation was then manually delineated and clipped. Finally, the excavated volume was computed as a 2.5D volume on a 1 mm resolution grid.

Automatic post-processing method

Accepted Article

After automatically trimming the borders of the excavations, dense 3D point clouds resulting from the flash lidar and SfM acquisitions were processed in several steps consisting of a combination of calls to CloudCompare 2.10.1 functions using R software and raster processes using the 'raster' R library. Excavated holes were detected using cloud-to-cloud distances between 3D point cloud before and after excavation, considering the minimal areas of the holes and the level of noises of the methods (Figure 2). The complete method is available as an R code in the open access repository Zenodo at the following URL (https://doi.org/10.5281/zenodo.4036423), and a sample dataset is available at the following URL (https://doi.org/10.5281/zenodo.4036313). <Figure 2 here>

Comparison of operating times

The long operating time of core sampling is mainly due to the high number of samples required to improve the statistical significance of the measured bulk density. The other methods require similar, lower operating times. However, the rubber balloon method in the field is also timeconsuming compared to the sensor-based methods (Table 4). <Table 4 here>

Statistical analyses

First, manual and automatic processing for volume estimation were compared two-by-two (SfM and flash Lidar) using a linear model. Second, taking the automatic volume processing as the basis for volume calculations, volumes issued from SfM method were compared to those issued from rubber balloon and flash lidar using a linear model and the classical figures of merit (R 2 , RMSE, intercept and slope) of the model as descriptors of the accuracy of the correlation between variables. Finally, all densities calculated from the automatic volume processing were compared using ANOVA considering the methods, horizons and replicates as explanatory variables. Reliability of the SfM method to the others was tested using differences in means and Tukey's HSD post hoc criterion.

Accepted Article

This article is protected by copyright. All rights reserved.

Results

Representativity of the dataset in terms of volumes and densities

Range of raw density ρ b

Based on the classical core sampling reference method, the studied horizons corresponded to a large observed range of ρ b [1.04 -1.71] mainly due to their various structural arrangements, organic carbon content and porosities. This dataset covered a large part of the bulk density measured in national databases (for example [START_REF] Hollis | Empirically-derived pedotransfert functions for predicting bulk density in European soils[END_REF]. The highest ρ b values measured with core sampling were in Mauguio_1_H2, due to tillage and wheeling degradations. The highest ρ b values (a mean of 2.01, measured in Mauguio_3_H2 with flash lidar and SFM) could not be measured with core sampling due to high pebble and gravel content (the highest : 73% in mass).

The lower ρ b values corresponded to the organic horizon Angles_H1 (1.04). The ρ bFE were lower than the ρ b values due to the high content of coarse fragments in some horizons.

Range of measured volumes

The distribution of the volumes measured with rubber balloon, SFM and flash lidar demonstrate a large range from 100 to 1400 cm 3 . The high number of measurements around the median (385 cm 3 ) corresponded to the rubber balloon excavations.

Comparison of the volume estimation from SfM, flash lidar, and rubber balloon

The automatic post-processing method used to estimate the volume from sensor-based methods was compared to the manual post-processing method. The automatically estimated volumes from both SfM and flash lidar were extremely similar to those estimated by the manual approach (R 2 =0.998, RMSE of 11 cm 3 ). Therefore, automatic post-processing was used in the rest of the analysis.

In the case of the investigation of the horizons from above, the excavations were sized for the rubber balloon method, and exactly the same holes were scanned by the flash lidar and SfM.

Hence, for these samples, exactly the same volumes from the three methods were measured and compared. The volumes measured by rubber balloon were systematically greater than the volumes measured by SfM (intercept of the linear model = 12 cm 3 , P<0.05 Student's t-test) (Figure 3a). Conversely, no bias was observed between the volumes from flash lidar and SfM (Figure 3b). However, the relationship between the volumes was noisy, especially at larger volumes, with an RMSE of 36 cm 3 .

Accepted Article

This article is protected by copyright. All rights reserved.

<Figure 3 here >

Determination of the bulk densities at the horizon scale <Table 5 here>

A comparison between the methods was conducted on the basis of the ρ bFE (figure 4 andtable 5).

The ρ bFE values measured from sensor-based methods were similar regardless of the soil characteristics, except for Lavalette_H2 (Tukey's HSD test, P<0.01). However, the standard deviation of the flash lidar data was systematically higher than those of SFM. The significant noisy measurement of the volume (an RMSE of 36 cm 3 ) directly corresponds to a high standard deviation of the resulting BD (more than 0.05 for a volume of 1000 cm 3 ). The comparison between classical methods and SFM showed differences in relation to the coarse fragments. For low content of coarse fragments (i.e. <15% in mass), The ρ bFE values measured from core sampling were significantly the same of those from SFM (Tukey's HSD test, P<0.01), except for Mauguio_1_H3 and Mauguio_2_H2. Moreover, the standard deviations of ρ bFE from both core sampling and SfM measured within each horizon were also similar. The ρ bFE values measured from rubber balloon was systematically lower than that measured from SFM and significantly lower for Lavalette_H2, Mauguio_3_H3, and Pech_2 (Tukey's HSD test, P<0.01). <figure 4 here> Conversely for high content of coarse fragments, the ρ bFE measured from core sampling were always lower than the other methods, regardless of the other soil properties (Tables 5 and2).

Finally, the variability of the measurements was high whatever the methods. The flash lidar resulted to the highest standard deviation measured in the Angles_H1 horizon and Mauguio_3_H2 horizon (the highest coarse fragment content).

Accepted Article

Discussion

The ten horizons were investigated both on the side of the pit (core sampling, SfM, flash LiDAR) and from above (rubber baloon, SfM, flash LiDAR). This methodology provided an original dataset which opens the discussion on the applicability and efficiency of each method. To the best of our knowledge, this study is the first experimental comparison of four different methods for soil bulk density measurement in a wide range of soil types and horizon depths.

Accuracy of SfM for soil bulk density measurement

The acquisition of the excavated volume is key for soil bulk density measurements.

Reference core sampling uses a fixed volume to extract the bulk soil. The other methods, except gamma rays, depend on an excavated volume of soil. Of the SfM, rubber balloon and flash lidar methods, SfM provided the greatest improvement over the core sampling method in terms of its measurements of bulk density. This result is mainly due to the high accuracy of the volume measurement, which had a very low standard deviation. In line with the results from [START_REF] Bauer | A photogrammetric method for calculating soil bulk density[END_REF] from the soil surface, our study, applied to soil horizons, helped to cover all components of pedological studies. In addition, the measure of the volume with the rubber balloon was systematically higher than that measured by the SfM, leading to lower resulting bulk densities. This was already observed in many studies (Ship and Matelski, 1965;[START_REF] Yoro | Les méthodes de mesure de la densité apparente. Analyse de la dispersion des résultats dans un horizon donné[END_REF][START_REF] Andraski | Balloon and core sampling for determining Bulk Density of alluvial desert soil[END_REF][START_REF] Page-Dumroese | Comparison of methods for determining bulk densities of rocky forest soils[END_REF]. The main proposed reason is that the size of the excavation is larger than the volume sampled with core sampling, and allows to better taking into account the structural porosity. In our study, exactly the same excavated volume was estimated with the rubber balloon, flash lidar and SfM methods. Therefore, our findings suggest that the volumes of the excavations may be systematically overestimated when measured with the rubber balloon, probably due to the overestimation of the first volume corresponding to the void between the edge of the device and the soil surface. Moreover, the pressure of the rubber balloon probably flattens out the edge of the excavation. These hypotheses are also corroborated by the bulk densities measured with core sampling, which showed no bias when compared with the bulk densities measured by SfM both from above and from the side of the pit. Our results may hence indicate that, in addition to the already known source of bias linked to larger investigated volumes, the rubber balloon method may have inherent flaws in the volume estimation itself.

Moreover, the volume measured with rubber balloons seems to be noisier than that measured by SfM within the same horizon, especially when the gravel and pebble contents increase. The

Accepted Article

balloon may not truly match the edges of the excavations, although the shape follows smooth edges during the excavation. The coarse gravels and pebbles may modify the contact between the balloon and the excavation edge and result in an increase in noise.

Excavation methods are generally used in the literature when the gravel and stone content of the soil exceeds 10% by mass. The different horizons tested in our study show the satisfactory applicability of the SfM method in these cases. Indeed, the total amount of gravel and stone varied from 0 to 73% (in mass), and the SfM method was not impacted by the gravel content. Similarly, the rubber balloon method is not recommended in the case of sandy textures due to the collapse of the excavation [START_REF] Andraski | Balloon and core sampling for determining Bulk Density of alluvial desert soil[END_REF]. In our study, the SfM method succeeded in the sandy Pech_H2 horizon (79% of sand). However, it must be noted that the bulk density was slightly overestimated compared to that estimated by core sampling. Kutilek and Nielsen (1994) discussed about the representative elementary volume (REV) in Hydrodynamic and calculated a REV minor than 100 cm 3 for poor structured soils and larger for well-structured soils. The core sampling with a fixed investigated volume is limiting in order to characterise structured soils. The size of the excavation for SfM measurement may be adapted to each situation to take into account the different specific porosities (small cracks, structural vughs, pebble size). This point was already made by [START_REF] Bauer | A photogrammetric method for calculating soil bulk density[END_REF] in the case of small soil surface cracks that could not be characterized by classical core sampling. [START_REF] Scanlan | Three-dimensional scanning for measurement of bulk density in gravelly soils[END_REF] noted the value of sizing the excavation to characterize different representative volumes. For example in Mauguio_1_H3 with a high porosity combined to a very clear structure, the volume investigated with the core sampling were lower than those from SFM and flash lidar (more than 1000 cm 3 ) and probably did not include the structural porosity. Consequently, the ρ bFE values from core sampling were higher than from the other methods. Finally, an additional advantage of SfM is its applicability regardless of the orientation of the soil face, along the side of the pit or from above.

Displaying SfM from the field to the lab

In addition to the demonstrated accuracy and precision of the volume estimation both from above and from the side of the pit, a great advantage of the proposed SfM method is the very low time needed in the field (Table 4), as already noted in previous studies [START_REF] Bauer | A photogrammetric method for calculating soil bulk density[END_REF][START_REF] Scanlan | Three-dimensional scanning for measurement of bulk density in gravelly soils[END_REF]. Moreover, the volume is directly measured in the field and any undisturbed samples need to be returned to the laboratory for fine volume measurement, as proposed SfM methods from [START_REF] Moret-Fernandez | Applicability of the photogrammetry technique to determine the volume and the bulk density of small soil aggregates[END_REF] or from [START_REF] Whiting | Rapid bulk density measurement using mobile device photogrammetry[END_REF]. Field work is often performed under time constraints. The representativity of the sampling with only 5 samples (a total of 1500 Accepted Article cm 3 on average) is at least comparable to that of the core sampling with an equivalent of 15 samples. However, the sampling must follow important recommendations mainly: (i) maintain a smooth shape of the excavation, (ii) be careful to maintain a good lighting of the image without potential moving objects , and (iii) prevent external soil particles from falling during the excavation.

Almost all image data processing is automated, with the exception of the application of the Time-SIFT method. Currently, operator intervention is still needed between the two main steps (image alignment and dense matching) of Time-SIFT processing to separate the so-called before and after datasets prior to the computing of the dense 3D point clouds. This was not a limiting factor in the case of this study, which suggested that this study had a manageable amount of samples (several dozens), but may become a true limitation for the normal operational use of SfM-based soil bulk density measurements. There is hence a strong interest in developing a Time-SIFT plugin for Photoscan so that the whole image processing process would be fully automated, from image dataset provision to individual volume estimations. This work is already ongoing but is out of the scope of this study. Other further work may assess the potential for the use of simple sensors such as those on smartphones for image acquisition and further 3D processing. Even if the lens and sensor quality of smartphones is lower than that of single-lens reflex digital cameras, the quality of the latest phone cameras and the availability of multi-lens smartphones may provide enough quality to allow for accurate and precise volume measurements. A test realised on one excavation (data not shown) results in an excellent agreement between camera-based and phone-based SfM.

The promise of flash lidar

The soil bulk densities measured from the flash lidar data were generally significantly different from the values estimated with other methods. A prominent finding is the large standard deviation of flash lidar data compared to that from the SfM method. This calls into question the accuracy of the flash lidar method, which also depends on a large set of local environmental conditions, such as wind and scene illumination. Compared with that in our results, [START_REF] Scanlan | Three-dimensional scanning for measurement of bulk density in gravelly soils[END_REF] achieved better accuracy on BD measurements while using a Kinect™ sensor. The main factor explaining the difference between our findings and the findings of [START_REF] Scanlan | Three-dimensional scanning for measurement of bulk density in gravelly soils[END_REF] is most likely the sensor resolution. The Kinect sensor used by these authors is a 2 M pixel sensor (1920x1080), while the Pico Flexx used in our study is only a 38 k pixel sensor (224 x 171). With such a small number of pixels, and especially for the case of the investigation on the side of the Accepted Article pit, sampling of the 3D surface may not be sufficient to correctly depict the geometries before and after excavation. On the other hand, [START_REF] Polyakov | An improved excavation method for measuring bulk density of rocky soil using terrestrial LiDAR[END_REF] used a professional-grade terrestrial laser scanner whose accuracy and precision specifications are totally different from the specifications of the lightweight flash lidar such as the ones used in the study of [START_REF] Scanlan | Three-dimensional scanning for measurement of bulk density in gravelly soils[END_REF] and in this study. However, their work confirms the interest in using laser-based sensors for volume measurements, and our three studies definitely demonstrate the strong potential of 3D sensors for soil studies. Considering the differences among the sensors used in our three studies, it is still hard to determine whether continuous flash lidar scanning (as proposed by [START_REF] Scanlan | Three-dimensional scanning for measurement of bulk density in gravelly soils[END_REF] is needed to obtain the required precision and accuracy or whether a single point of view with time averaging (as we propose in our study), but with a better sensor would be sufficient.

Regardless, the main advantages of using such lightweight range imaging sensors are their easy handling without prior specific knowledge and the large frame (more than 2 m²) recorded during just a few seconds. These devices also connect directly to a simple smartphone and do not need a single-lens reflex camera, which is advised for the SfM. Finally, it is expected that better lightweight flash lidar devices will become available and would allow the use of this method to obtain the required accuracy and precision for the volume estimations.

The bulk density of fine earth from the stony soils

The studied dataset showed large differences with the measurement of bulk density in relation to the coarse fragment contents. Core sampling was considered as a reference for the measurement of bulk density in numerous studies. In case of low coarse fragment contents (i.e. <15%), the sampling was possible with the core sampling and a simple correction is sufficient to study the bulk density of the fine earth. Restinclieres_H2 and Angles_H1, with respectively 24 and 36 % in mass of coarse fragment contents, presented lower bulk densities of the fine earth than expected with the core sampling method. During the sampling in the field, the fine gravels and gravels were easily included in the core stainless steel and did not significantly deflect the driving of the core.

However, the stones modified the driving and the structure of the sample. The resulting sample consequently contained lower quantity of fine earth. Conversely, the excavation methods undisturbed the structure of the fine earth whatever the gravel and stone content. The comparison with the core sampling was biased and conducted to an underestimation of ρ bFE with the core sampling. In case of very high stone content (Mauguio_3_H2), the core sampling was

Accepted Article

even not applicable at all. The sensor-based excavation methods denoted low and noisy ρ bFE values, which probably depend on the irregular voids between the stones.

Stony soils (i.e. soils with more than 30 % in mass of coarse fragments) represent a large part of the Mediterranean soils and more generally of the European countries (The soil map of Europe, 2014). New studies focusing on the relation between stony soils and water retention show the importance of stones and the relation between stones and fine earth (for example, [START_REF] Korboulewsky | Plants use water in the pores of rock fragments during drought[END_REF]). The structural arrangements within the fine earth in stony soils needs to be studied and the new proposed excavation methods may improve future datasets.

Conclusion

Although BD is a key property in soil science, the classical methods for measuring BD in the field are time-consuming and present some limitations, especially in stony soils. The aim of the present work was to assess SfM as an all-terrain method capable of overcoming the limitations of conventional methods for BD measurements, regardless of the stone content. The SfM method provides accurate measurements with significantly less time in the field than the traditional method, regardless of the soil conditions. Moreover, the representativeness of the measurement is larger than that of the reference core sampling and is adjustable according to the size of the different soil porosities. The volume of excavations measured from rubber balloons, which are often used in the case of stony soils, is significantly higher than the volume measured from the SfM. This result confirms the bias of the rubber balloon method, which is often discussed in the literature. At this stage, lightweight flash lidar is not sufficient to measure BD with high accuracy.

However, the flash lidar method is very accessible, and technical advances may provide more accurate devices in the future. Finally, the authors recommend the use of SfM to measure BD.

Acknowledgements

This study was directly supported by the LISAH research team. The authors thank N. Saurin for his kind help at the Pech Rouge site, the DIASCOPE team and Y. Blanca for their appreciated help at the Mauguio site. *This time is lower mainly due to the ease of vertical sampling. In the case of horizontal sampling (rubber balloon), the excavation is time-consuming. Horizons were ranked by their coarse fragment content, from low (Lavalette_H4) to high (Mauguio_3_H2).

Accepted Article

Accepted Article

Bold cells with * represent significantly equal means for each line (Tukey's HSD based on the Student distribution, at P<0.01) comparing to SFM.

Table 1 :

 1 Morphological properties of the investigated soil horizons

	Investigated soil	WRB	minimum	maximum	Colour	Structure	porosity
	horizon		depth	depth
		Soil type
			(m)	(m)
						Type*	Size (mm)	pore level	type
	Lavalette_H2	Calcisol	0.18	0.45	2.5Y44	SB/MA	20	medium/low	vughs
	Lavalette_H4		0.60	1.90	2.5Y54	SB	30	very high	channels
	Mauguio_1_H2	Fluvisol	0.65	1.00	2.5Y43	PR/SB	100/30	medium	vughs
	Mauguio_1_H3		1.00	1.40	2.5Y64	SB	30	high	vughs
	Mauguio_3_H2	Rhodic	0.45	0.75	5YR56	AB	10	low	-
		luvisol
	Mauguio_3_H3		0.75	1.20	5YR44	AB	40	medium	-
	Mauguio_2_H2	Fluvisol	0.40	0.80	10YR54	SB/MA	30	low	-
	Pech_H2	Sandy	0.35	0.60	10YR44	SB	30	low	-
		Cambisol
	Angles_H1	Cambic	0	0.35	10YR22	GR	10	medium	vughs
		Umbrisol
	Restincliere_H2	Calcisol	0.10	0.25	10YR54	SB	20	medium	channels
	*according to FAO Guidelines for soil description (2006) : SB:subangular blocky; AB:angular blocky; MA:massive; PR:prismatic;
	GR:granular.

Table 2 :

 2 Main physico-chemical properties of the investigated soil horizons

	Accepted Article	Investigated horizon Lavalette_H2 Lavalette_H4 Mauguio_1_H2 Mauguio_1_H3 Mauguio_3_H2 Mauguio_3_H3 Mauguio_2_H2 Pech_H2 Angles_H1 Restincliere_H2 soil	fine gravel (%mass) 0 0 1 1 11 3 2 4 14 8	gravel/stones (%mass) 2 0 2 0 62 3 3 8 22 16	density of the coarse fragments 2.5 2.5 2.5 2.5 2.3 2.3 2.5 2.5 1.9 2.6	A (%) 24 27 26 23 48 45 19 10 12 20	L (%) 47 47 42 39 27 33 32 11 24 42	S (%) 29 26 32 38 25 22 49 79 64 38	Organic carbon (g/kg) 6.2 4.1 6 4.5 3.5 2.4 3.7 5 41 6.7	Calcium carbonate (g/kg) 579 734 355 367 0 0 348 306 0 749

This article is protected by copyright. All rights reserved.

Table 4 :

 4 Average human operating time (in min) of each method for one horizon with sufficient replications

	Time (min)

Table 5 :

 5 Mean soil bulk densities of fine earth (ρ bFE ) measured by the different methods within each

	investigated horizon
	horizons	Coarse	Sensor-based methods	Classical methods
		fragments
		(% in mass)	SfM	Flash lidar	Rubber balloon	Core sampling
	Lavalette_H4	0	1.30 ± 0.06 (10)*	1.28 ± 0.14 (10)*	1.22 ± 0.07 (5)*	1.33 ± 0.06 (15)*
	Mauguio_1_H3	1	1.53 ± 0.04 (10)*	1.45 ± 0.10 (10)*	1.47 ± 0.05 (5)*	1.59 ± 0.03 (15)
	Lavalette_H2	2	1.63 ± 0.04 (10)*	1.57 ± 0.04 (10)	1.55 ± 0.04 (5)	1.65 ± 0.04 (15)*
	Mauguio_1_H2	3	1.67 ± 0.04 (10)*	1.62 ± 0.10 (8)*	1.65 ± 0.09 (5)*	1.70 ± 0.07 (4)*
	Mauguio_2_H2	5	1.64 ± 0.03 (10)*	1.56 ± 0.09 (10)*	1.57 ± 0.11 (5)*	1.69 ± 0.05 (15)
	Mauguio_3_H3	6	1.59 ± 0.03 (10)*	1.68 ± 0.18 (10)*	1.51 ± 0.06 (5)	1.58 ± 0.04 (15)*
	Pech_H2	12	1.66 ± 0.05 (4)*	1.70 ± 0.03 (2)*	1.46 ± 0.08 (5)	1.

61 ± 0.07 (15)* Restinclieres_H2 24 1.43 ± 0.09 (10)* 1.41 ± 0.13 (10)* 1.39 ± 0.13 (5)*

						1.28 ± 0.13 (15)
	Angles_H1	36	1.03 ± 0.07 (10)*	1.14 ± 0.23 (10)*	0.93 ± 0.08 (5)*	0.88 ± 0.05 (14)
	Mauguio_3_H2	73	1.45 ± 0.22 (8)*	1.99 ± 1.42 (8)*	1.12 ± 0.18 (4)*	-
	Mean ± standard deviation (number of samples)

This article is protected by copyright. All rights reserved.

Data availability statement:

R code https://doi.org/10.5281/zenodo.4036423

Conflict of interest: None

Accepted Article

This article is protected by copyright. All rights reserved.