Mechanical diffusion in grease ice stirred by gravity waves

TL;DR

This paper investigates how gravity waves can induce mechanical diffusion in grease ice through hydrodynamic interactions, modeling the ice as a deformable porous medium and analyzing the resulting effective diffusivity.

Contribution

It introduces a two-fluid model for grease ice under gravity waves, highlighting the potential for enhanced diffusion due to mechanical interactions among ice crystals.

Findings

Effective diffusivity can surpass contact-induced dislocation values.

Hydrodynamic interactions lead to interstitial flows similar to diffusion in fixed beds.

The model predicts parameter ranges for significant diffusion enhancement.

Abstract

The possibility of hydrodynamic diffusion in a model of grease ice stirred by the velocity field of a gravity wave is explored. It is argued that mechanical interactions among ice crystals can induce disturbances in the fluid velocity - in the form of interstitial flows - analogous to those leading to diffusion in fixed beds. A two-fluid description of the system is introduced, in which the ice matrix is treated as a deformable porous medium. Depending on the range of parameters, the effective diffusivity of the medium can exceed the value which would be obtained by only keeping into account particle dislocation induced by contact interactions.

Equations

$${\kappa }^{wind}\approx (u_{\ast }L{)}^2/{\nu }_{grease}\approx 10^{-8}{\mathrm{m}}^2/s.$$

$${\kappa }^{wave}\approx (ϵL{)}^2\omega \approx 10^{-8}{\mathrm{m}}^2/s.$$

$$(1-\bar{C})[{\rho }_l\frac{\partial {\mathbf{U}}_l}{\partial t}+\nabla P-\mu {\nabla }^2{\mathbf{U}}_l]=\Gamma ({\mathbf{U}}_s-{\mathbf{U}}_l)+(1-C){\rho }_l\mathbf{g},$$

$$\bar{C}[{\rho }_s\frac{\partial {\mathbf{U}}_s}{\partial t}+\nabla P-\mu {\nabla }^2{\mathbf{U}}_l]=-\nabla P_s+{\mu }_s^B\nabla \nabla \cdot {\mathbf{U}}_s+{\mu }_s{\nabla }^2{\mathbf{U}}_s-\Gamma ({\mathbf{U}}_s-{\mathbf{U}}_l)+C{\rho }_s\mathbf{g},$$

$$\frac{\partial \tilde{C}}{\partial t}+\bar{C}\nabla \cdot {\mathbf{U}}_s=0,$$

$$\nabla \cdot [(1-\bar{C}){\mathbf{U}}_l+\bar{C}{\mathbf{U}}_s]=0,$$

$${\bar{P}}_s(z)={\bar{P}}_s(0)+\bar{C}g({\rho }_l-{\rho }_s)z\mathrm{and}\bar{P}(z)=\bar{P}(0)-g{\rho }_{l}z.$$

$$(1-\bar{C})[{\rho }_l\frac{\partial {\mathbf{U}}_l}{\partial t}+\nabla \tilde{P}-\mu {\nabla }^2{\mathbf{U}}_l]=\Gamma ({\mathbf{U}}_s-{\mathbf{U}}_l),$$

$$\bar{C}[{\rho }_s\frac{\partial {\mathbf{U}}_s}{\partial t}+\nabla \tilde{P}-\mu {\nabla }^2{\mathbf{U}}_l]=-\beta \nabla \tilde{C}+{\mu }_s^B\nabla \nabla \cdot {\mathbf{U}}_s+{\mu }_s{\nabla }^2{\mathbf{U}}_s$$

$$-\Gamma ({\mathbf{U}}_s-{\mathbf{U}}_l)+\tilde{C}({\rho }_s-{\rho }_l)\mathbf{g},$$

$$\beta =(\partial P_s/\partial C{)}_{C=\bar{C}}$$

$$\tilde{C}(\mathbf{r},t)=\tilde{C}(z){\mathrm{e}}^{\mathrm{i}(kx-\omega t)},\tilde{C}(z)={\tilde{c}}_p{\mathrm{e}}^{kz}+\sum _n{\tilde{c}}_n{\mathrm{e}}^{{\alpha }_{n}z},\dots ,$$

$${\mathbf{U}}_s=-\nabla \Phi -{\nabla }_{\perp }A,{\mathbf{U}}_l={\mathbf{U}}_s-\nabla \hat{\Phi }-{\nabla }_{\perp }\hat{A},$$

$${\tilde{c}}_p=0,{\tilde{c}}_n=-\frac{({\alpha }_n^2-k^2)\bar{C}}{\mathrm{i}\omega }{\varphi }_n\mathrm{and}{\varphi }_n=-(1-\bar{C}){\hat{\varphi }}_n.$$

$$(1-\bar{C})[\mathrm{i}\omega {\rho }_l+\mu ({\alpha }_n^2-k^2)](a_n+{\hat{a}}_n)-\Gamma {\hat{a}}_n=0,$$

$$\bar{C}[\mathrm{i}\omega {\rho }_s{a}_n+\mu (a_n+{\hat{a}}_n)]+{\mu }_s({\alpha }_n^2-k^2)a_n+\Gamma {\hat{a}}_n-\frac{g({\rho }_l-{\rho }_s)\bar{C}k}{\omega }(1-\bar{C}){\hat{\varphi }}_n=0,$$

$$\{[\mathrm{i}\omega {\rho }_l+\mu ({\alpha }_n^2-k^2)]\bar{C}(1-\bar{C})-\Gamma \}{\hat{\varphi }}_n+(1-\bar{C})p_n=0,$$

$$\{[-\mathrm{i}\omega {\rho }_s\bar{C}-({\mu }_s^B+{\mu }_s-\bar{C}\mu -\frac{\beta \bar{C}}{\mathrm{i}\omega })({\alpha }_n^2-k^2)$$

$$-\frac{g({\rho }_l-{\rho }_s)\bar{C}{\alpha }_n}{\mathrm{i}\omega }](1-\bar{C})+\Gamma \}{\hat{\varphi }}_n+\bar{C}{\tilde{p}}_n=0.$$

$$(1-\bar{C})[\mathrm{i}\omega {\rho }_l({\varphi }_p+{\hat{\varphi }}_p)+p_p+{\rho }_{l}g{\eta }_l]-\Gamma {\hat{\varphi }}_p=0,$$

$$\bar{C}[\mathrm{i}\omega {\rho }_s{\varphi }_p+p_p+{\rho }_{s}g{\eta }_s]+\Gamma {\hat{\varphi }}_p=0,$$

$$-\mathrm{i}\omega {\eta }^{s,l}=U_z^{s,l}{|}_{z=0}.$$

$$\gamma =1+\frac{{\mu }_s^B}{{\mu }_s},\rho =\frac{{\rho }_s}{{\rho }_l},\mathrm{and}{\nu }_g=\frac{k_{\omega }^{3/2}{\mu }_s}{{\rho }_l{g}^{1/2}},$$

$$\mu ,{\nu }_g,\bar{C}\ll 1.$$

$$\mathrm{i}{a}_n-\Gamma {\hat{a}}_n\simeq 0,$$

$$[\mathrm{i}\rho \bar{C}+{\alpha }_n^2]a_n+\Gamma {\hat{a}}_n-(1-\rho )\bar{C}{\hat{\varphi }}_n\simeq 0,$$

$$\Gamma {\hat{\varphi }}_n-p_n\simeq 0,$$

$$[\mathrm{i}(1-\rho ){\nu }_g^{-1/2}\bar{C}{\alpha }_n-(\gamma +\mathrm{i}\beta ){\alpha }_n^2+\Gamma ]{\hat{\varphi }}_n+\bar{C}{p}_n\simeq 0$$

$${\alpha }_i\simeq \sqrt{-\mathrm{i}},{\hat{a}}_i\simeq \mathrm{i}{a}_i/\Gamma ,{\hat{\varphi }}_i\simeq p_i\simeq 0.$$

$${\alpha }_c\simeq \frac{\mathrm{i}(1-\rho )\bar{C}\pm \sqrt{4(\gamma +\mathrm{i}\beta )\Gamma {\nu }_g-(1-\rho {)}^2{\bar{C}}^2}}{2(\gamma +\mathrm{i}\beta )\sqrt{{\nu }_g}},$$

$$p_c\simeq \Gamma {\hat{\varphi }}_c,a_c\simeq \frac{(1-\rho )\bar{C}}{{\alpha }_c^2}{\hat{\varphi }}_c,{\hat{a}}_c\simeq \frac{\mathrm{i}(1-\rho )\bar{C}}{{\alpha }_c^2\Gamma }{\hat{\varphi }}_c.$$

$$[U_s-U_l{]}_{z=0}\simeq {\hat{\varphi }}_p+{\alpha }_c{\nu }_g^{-1/2}{\hat{\varphi }}_c=0,$$

$${\eta }_s\simeq {\eta }_l\simeq -\mathrm{i}k{\varphi }_p-a_i+\mathrm{i}{\alpha }_c{\nu }_g^{-1/2}{\hat{\varphi }}_c.$$

$$\mathrm{i}(1-k){\varphi }_p+{\tilde{p}}_p-a_i+\Gamma {\nu }_g^{-1/2}{\alpha }_c{\hat{\varphi }}_c=0,$$

$$\bar{C}[\mathrm{i}\rho (1-k){\varphi }_p+{\tilde{p}}_p-\rho a_i]-\Gamma {\nu }_g^{-1/2}{\alpha }_c{\hat{\varphi }}_c=0,$$

$${\hat{\varphi }}_c\simeq \frac{(1-\rho )\bar{C}{\nu }_g^{1/2}}{\Gamma {\alpha }_c}{\tilde{p}}_p,{\tilde{p}}_p\simeq -\mathrm{i}(1-k){\varphi }_p+a_i.$$

$${\tau }_{zz}{|}_{z=0}\simeq -{\tilde{p}}_p+2\mathrm{i}{\nu }_g^{1/2}{\alpha }_i{a}_i+[(\mathrm{i}\beta +\gamma +1){\alpha }_c^2+2\mathrm{i}(1-\rho ){\nu }_g^{1/2}\bar{C}{\alpha }_c^{-1}]{\hat{\varphi }}_c=0,$$

Full text

1 Introduction

Flows in suspensions are often characterized by transport phenomena that originate from the interaction of the particles among themselves and with the flow. In dilute suspensions, the particle interactions are mediated by hydrodynamic forces. In more concentrated suspensions, the contribution from direct contact of the particles is dominant. Effective diffusivities greatly exceeding the values expected from molecular effects are often observed, even in the absence of turbulence [1, 2].

A medium of particular interest in the geophysical context is grease ice, that is the dense slurry of ice crystals often observed at the water surface in polar seas during periods of strong wind and freezing temperatures [3]. Consolidation of the slurry to pack ice rests on efficient removal of the heat and salt produced in freezing. Due to the high viscosity of the medium, turbulence is severely hindered. Mechanical diffusion induced by stirring of the ice by the wind and by gravity waves, on the contrary, may play an important role.

A possible mechanism is shear-induced diffusion [4, 5, 6]. It is possible to make some estimate in the two cases diffusion is generated by wind stirring and by gravity waves.

The wind will generate in the ice layer a shear of strength $a_{s}\approx(u_{*})^{2}/\nu_{grease}$, where $u_{*}$ is the friction velocity and $\nu_{grease}$ is the effective viscosity of grease ice. In a concentrated suspension, such as grease ice, particle diffusion will arise from collisions, rather than from hydrodynamic interaction between particles. Collisions will occur between ice crystals travelling along flow lines separated by a typical crystal size $L$, corresponding to a collision cross section $\sigma\approx L^{2}$, and to a typical relative velocity in collision $u_{L}\approx La_{s}$. In concentrated conditions, the number density of the crystals is $n\approx L^{-3}$ and the particle displacement in a collision is $\approx L$. Hence, the collision rate will be $a_{coll}\approx n\sigma_{coll}u_{L}\approx a_{s}$, and the wind contribution to the diffusivity $\kappa_{wind}\approx a_{s}L^{2}$. Taking a wind speed of about 10 m/s, corresponding to a friction velocity $u_{*}\approx 0.01$ m/s [7], and taking for the grease ice, $L\approx 1$ mm and $\nu_{grease}\approx 0.01\ {\rm m^{2}/s}$ [3], would give

[TABLE]

The contribution to diffusion from gravity waves can be estimated in a similar way. A crystal in the field of a gravity wave will be able to “prod” a nearby crystal only during the compressible phase of the wave at the crystal location. The resulting displacement will be $S\approx u_{L}/\omega\approx La_{s}/\omega$, where $\omega$ is the wave frequency, $a_{s}\approx\epsilon\omega$ is the strain intensity, $\epsilon=kU/\omega$ is the dimensionless wave amplitude, and $k$ and $U$ are the wavenumber and the characteristic fluid velocity of the wave. Hence $S\approx\epsilon L$. In concentrated conditions, a crystal will make collision with a neighbour about once in a period, and if the degree of irreversibility in collision is sufficient [8], a contribution to the diffusivity, $\kappa^{wave}\approx S^{2}\omega$, will be generated. Taking typical values for ocean waves, $\omega\approx 1\ {\rm s^{-1}}$ and $\epsilon\approx 0.1$, would give

[TABLE]

If molecular diffusion is small, tracers will move with the liquid around the particles, and have the same diffusion properties of the particles, described by Eqs. (1) and (2). The diffusivities in Eqs. (1) and (2) can then be compared with the heat and salinity diffusivities in salt water: $\kappa_{T}\simeq 1.4\times 10^{-7}\ {\rm m^{2}/s}$ and $\kappa_{S}\simeq 7.4\times 10^{-10}\ {\rm m^{2}/s}$. The mechanical contribution to diffusion appears to be relevant only for salinity.

The estimates in Eqs. (1) and (2) disregard the possible impact of relative motions between phases, which may lead to diffusion phenomena analogous to those generated in flows in fixed particle beds [9, 10]. Relative motion of the phases is going to be generated in concentrated suspensions by the interplay between non-hydrodynamic and hydrodynamic (drag) forces on the particles. In creeping flow conditions, exact balance between mechanical and drag forces must exist, which corresponds, macroscopically, to balance between the drag forces by the interstitial flow and the resistance to deformation in the solid matrix. In diluted suspensions, no mechanical interparticle interactions obviously exist, and balance of solid and fluid stresses occurs at the particle scale, that is also the scale of the fluid velocity disturbances.

In grease ice, non-hydrodynatic stresses are generated dynamically by the continuous formation and breaking of bonds between ice crystals induced by the flow, and allow the medium to persist in a fluid state.

Note that in order for a crossflow to be present, a pressure gradient must exist in the liquid phase. Such pressure gradients are absent in the homogeneous shear layer generated by the wind stress. They are generated dynamically, instead, in the field of a gravity wave. The question remains whether the contribution to diffusion is significant or not. Answering this question is the goal of the present paper.

To study the problem, a two-phase approach is adopted. No attempt is made to determine the macroscopic stresses of the network from microscopic behaviors, rather, an isotropic viscous response to deformation is assumed, with viscosity extrapolated from that of real grease ice. Knowing the relative velocity of the phases allows to obtain diffusivity estimates based on the characteristic size of the particles in suspension. Not more than an order of magnitude estimate is provided: grease ice is a complicated mixture of mm and sub-mm crystals in the form of platelets, needles and other shapes, with a solid volume fraction $\bar{C}\approx 0.1$ [7, 11]. Analytical treatment of the problem is therefore very difficult. More precise results on diffusion in the field of a gravity wave are obtained in the Appendix, in the idealized case of a dilute fiber network.

The analysis predicts significant diffusivity enhancement with respect to the estimates in Eqs. (1) and (2), with an amplification factor proportional to the ratio between the distance traveled by a tracer in the solid matrix and the crystal scale. The key condition for amplification is that the effective viscosity of the slurry is large compared to that of the interstitial liquid. The diffusivity dependence on the solid volume fraction is weaker.

The paper is organized as follows. In Sec. 2, the macroscopic equations describing the two-phase system are derived. In Sec. 3, the equations are applied to the calculation of the relative motion between phases. In Sec. 4, the contribution to diffusivity is evaluated. Section 5 is devoted to conclusions. Analytical results for diffusion in a dilute random fiber bed are presented in the Appendix.

2 Two-phase model

The dynamics of gravity waves in grease ice covered ocean is studied by modelling the ocean as an infinitely deep uniform slurry. As the thickness of grease ice in the ocean, and the depth of the region where wave dissipation takes place, are usually comparable [12, 13], the model is going to be rough, but not completely inaccurate.

A macroscopic two-phase description is adopted, with the solid and liquid phases obeying separate mass and momentum conservation equations. No hypothesis is made for the moment on the microscopic structure of the medium, except for the fact that non-hydrodynamic interactions dominate the stress in the ice matrix.

To fix the ideas, take the $x$ axis along the direction of propagation of the wave, and the $z$ axis pointing upwards, with $z=0$ at the surface.

The momentum and continuity equations obeyed by the two phases read:

[TABLE]

where linearity has been imposed, tilde and overbar indicate fluctuating and mean quantities, $g\simeq 9.8\ {\rm m/s^{2}}$, $\Gamma$ and $C$ are the gravitational acceleration, the drag coefficient and the solid volume fraction, and $\rho_{n}$, ${\bf U}_{n}$ and $\mu_{n}$ are the mass density, the velocity and the dynamic viscosity of the liquid ($n=l$) and solid ($n=s$) phase. Similar equations have been used in [14] to model the flow in deformable porous media.

The terms $\nabla P_{s}$, $\mu_{s}^{\scriptscriptstyle{B}}\nabla\nabla\cdot{\bf U}_{s}$ and $\mu_{s}\nabla^{2}{\bf U}_{s}$ in Eq. (4) give the non-hydrodynamic part of the pressure and viscous forces. All three contributions arise from coarse graining of the particle mechanical interactions. Note that when these terms are zero, and the particles are neutrally buoyant, Eqs. (3) and (4) have solution ${\bf U}_{l}={\bf U}_{s}$. The mechanical pressure $P_{s}$ and the bulk viscosity $\mu_{s}^{\scriptscriptstyle{B}}$ account for the resistance to compression of the solid matrix. An isotropic expression is utilized for the bulk viscosity term even though such normal stresses are typically anisotropic in interacting particle suspensions [15]. As the flow of the individual phases is almost incompressible (see below), the expected error is small.

The terms $\nabla P$ and $\mu\nabla^{2}{\bf U}_{l}$ in Eqs. (3) and (4) give the hydrodynamic part of the pressure and viscous forces, with $\mu\nabla^{2}{\bf U}_{l}$ receiving contribution from the stress both in the liquid and in the particles (the Einstein correction in a dilute suspension). The hydrodynamic pressure $P$ is determined by the incompressibility condition Eq. (6). The mechanical pressure is taken to obey a constitutive relation in the form $P_{s}=P_{s}(C;\nabla{\bf U}_{s})$, and to rapidly vanish for $C>\bar{C}$, to allow hydrostatic equilibrium at $C\approx\bar{C}$.

Hydrostatic equilibrium for both phases is obtained by setting ${\bf U}_{l,s}=0$ in Eqs. (3) and (4). The result is, assuming slow variation of $\bar{C}$ throughout the domain,

[TABLE]

(Note that in the absence of a mechanical pressure, $\bar{P}(z)=\bar{P}(0)-g[\rho_{l}+\bar{C}(\rho_{s}-\rho_{l})]z$, buoyancy forces on the two phases, $-\bar{f}_{l}=\bar{f}_{s}=-\bar{C}(\partial_{z}\bar{P}+\rho_{s}g)=g(1-\bar{C})\bar{C}(\rho_{l}-\rho_{s})$, would be present, and would need to be balanced by the drag forces generated by settling of the particles to the surface.)

Subtraction of the hydrostatic component accounted for by Eq. (7), from Eqs. (3) and (4), gives

[TABLE]

where

[TABLE]

is the compressibility of the solid phase.

In the case of monocromatic waves, one can write for the fluctuating quantities:

[TABLE]

where subscripts $p$ and $n$ identify potential and non-potential components. The exponents $\alpha_{n}$ give the depth of the boundary layers induced by the two-phase dynamics at the water surface.

The velocities can be expressed in terms of velocity potentials [16]:

[TABLE]

where $\nabla_{\perp}=(\partial_{z},-\partial_{x})$. The fields $A$ and $\tilde{A}$ are chosen to contribute only to the vortical part of ${\bf U}_{l,s}$, and therefore do not have potential components $A_{p}=a_{p}{\rm e}^{kz}$ and $\hat{A}_{p}=\hat{a}_{p}{\rm e}^{kz}$. Such components could be eliminated anyway by a gauge transformation $\Phi\to\Phi-{\rm i}A_{p}$, $\hat{\Phi}\to\hat{\Phi}-{\rm i}\hat{A}_{p}$.

Substitution of Eq. (12) into Eqs. (5) and (6) gives

[TABLE]

Equations for the non-potential components $\tilde{p}_{n}$, $\hat{\phi}_{n}$, $a_{n}$ and $\hat{a}_{n}$ are obtained by taking $\nabla_{\perp}$ and $\nabla$ of Eqs. (8) and (9), and then using Eq. (13):

[TABLE]

The potential terms are obtained by integrating Eqs. (3) and (4) over $z$, and then exploiting Eqs. (7) and (13-17). The result is

[TABLE]

where the vertical displacement at the surface, $\eta^{s,l}$, obeys

[TABLE]

For weakly damped waves (which is the case in most situations of interest [17]), $k\simeq k_{\omega}=\omega^{2}/g$, and the dominant contribution to ${\bf U}_{l}\simeq{\bf U}_{s}$ is produced by $\phi_{p}$. The velocity scale of the waves is therefore $U=k_{\omega}|\phi_{p}|$.

3 Determining the crossflow

It is convenient to shift to dimensionless variables. The various quantities appearing in Eqs. (3-6) and below are rescaled as follows: $\omega t\to t$, $k/k_{\omega}\to k$, $\mu/\mu_{s}\to\mu$, $\sqrt{\mu_{s}/(\omega\rho_{l})}\alpha\to\alpha$, $\beta\bar{C}/(\omega\mu_{s})\to\beta$, $\Gamma/(\omega\rho_{l})\to\Gamma$, $pk_{\omega}/(\rho_{l}g)\to p$, $Uk_{\omega}/\omega\to U$. The velocity scale of the waves becomes $U=\epsilon$.

The following dimensionless parameters are introduced:

[TABLE]

and the following smallness conditions are imposed

[TABLE]

The stresses in the solid matrix are assumed to be of comparable magnitude; therefore $\gamma,\beta=O(1)$.

The conditions in Eq. (22) correspond, in the order, to the requirements:

• •

high viscosity of the medium compared to that of the liquid phase;

• •

waves that are only weakly damped [18, 13], in such a way that $|k-1|\ll 1$ and all velocity potentials are small except $\phi_{p}$;

• •

diluteness of the medium.

The first two conditions are easily satisfied by grease ice in the ocean. Assuming that all of the viscosity renormalization comes from mechanical interactions, in such a way that $\mu\approx\mu_{l}=\rho_{l}\nu_{l}$, $\nu_{l}\simeq 1.5\times 10^{-6}\ {\rm m^{2}/s}$, and taking for the viscosity in the solid matrix $\nu_{g}\approx\nu_{grease}\approx 0.01\ {\rm m^{2}/s}$ [3], would give $\mu\approx 10^{-4}$. In similar fashion, considering wavelengths in the range of $10-100$ m, would give for the dimensionless viscosity, $1.5\times 10^{-3}>\nu_{g}>5\times 10^{-5}$.

The last condition $\bar{C}\ll 1$ is only marginally satisfied, and is adopted for simplicity. An upper bound for the grease ice concentration is fixed by the transition from grease ice to compact ice, which takes place at $\bar{C}\approx 0.3$ [19].

Working with dimensionless variables and expanding to lowest order in the small quantities, Eqs. (14-17) take the form

[TABLE]

From Eqs. (23-26) an incompressible mode for the individual phases is identified: 111Note that in rescaled variables a factor $\nu^{-1/2}_{g}$ is extracted out of $\alpha_{n}$, in such a way that the $z$-derivative of non-potential fields $Q_{n}$ is $\partial_{z}Q_{n}=\nu_{g}^{-1/2}\alpha_{n}Q_{n}$.

[TABLE]

It is easy to recognize in this mode the standard behavior of the vortical component of the velocity field of a wave propagating in a viscous medium of viscosity $\nu_{g}$ [16, 13].

Equations (23-26) have an additional compressible mode with eigenmodes

[TABLE]

where only the root with a positive real part is admissible in an infinitely deep domain. The root’s inverse gives the thickness of the compressible boundary layer. From Eqs. (27) and (28), $|\alpha_{c}/\alpha_{i}|\sim\max(\Gamma^{1/2},\bar{C}\nu_{g}^{-1/2})\gg 1$, which implies that the compressible boundary layer is exceedingly thin.

From Eqs. (23-26) one can then write for the compressible component of the velocity potentials and of the pressure:

[TABLE]

To derive a dispersion relation, it is necessary to impose boundary conditions at the fluid free surface, and exploit Eqs. (18-20) to express the amplitude of non-potential modes in terms of potential ones. There are in total five unknowns: $k$, $p_{p}$, $\hat{\phi}_{p}$, $\hat{\phi}_{c}$ and $a_{i}$; $\phi_{p}$ remains arbitrary because of linearity. Thus, three boundary conditions are required. Two are well known from the study of gravity waves in viscous media [16], and are the requirement that the normal stress $\tau_{zz}\simeq-\tilde{P}-\tilde{P}_{s}+2\mu_{s}\partial_{z}U_{sz}+\mu_{s}^{\scriptscriptstyle{B}}\nabla\cdot{\bf U}_{s}$ and the tangential stress $\tau_{xz}\simeq\mu_{s}(\partial_{x}U_{sz}+\partial_{z}U_{sx})$ at the water surface are identically zero. A reasonable third boundary condition is that the interfaces with the atmosphere of the two phases move together, $[U_{s}-U_{l}]_{z=0}=0$, so that there are no sloshing phenomena at the water surface. Working with velocity potentials,

[TABLE]

and therefore,

[TABLE]

Substitution of Eqs. (30) and (31) into the equations for the potential terms, Eqs. (18) and (19), gives

[TABLE]

which allows to determine the amplitude of the compressible mode as a function of potential modes:

[TABLE]

The compressible mode is driven by buoyancy. The system of equations is closed by imposing the conditions on the normal and tangential stress at the surface:

[TABLE]

Substitution of Eq. (34) into Eqs. (35) and (36) gives finally

[TABLE]

One recovers the dispersion relation and the equation for the boundary layer structure of gravity waves in a homogeneous medium of kinematic viscosity $\nu_{g}\simeq\mu_{s}/\rho_{l}$ [16]. The suspension inherits the viscosity of the solid phase, as expected. Similar mechanisms have been invoked in [20] to explain the rheology of concentrated sphere suspensions. The compressible mode does not contribute to lowest order to the dynamics, as can be verified by back substitution of Eqs. (37) and (28) into Eqs. (35) and (36). This suggests that use of a boundary condition on $\eta^{s,l}$ could be avoided by simply setting $\hat{\phi}_{c}=0$, thus neglecting all contributions from the compressible mode. The thin compressible boundary layer at the surface would be replaced in this way by a finite gap between phases, whose amplitude can be evaluated from Eqs. (30) and (37): $\eta_{s}-\eta_{l}\simeq 2(1-\rho)\bar{C}\nu_{g}\Gamma^{-1}\phi_{p}$.

It is now possible to determine the relative motion of the phases. Inspection of Eqs. (27), (29), (30) and (37) allows to conclude that the dominant contribution is produced by the vortical incompressible mode $\hat{a}_{i}$. This leads to the final result

[TABLE]

The relative motion between phases is predominantly along the horizontal, and is confined in a region of thickness $\lambda_{\alpha}\approx\nu_{g}^{1/2}$ near the surface. The incompressible mode dominates the dynamics both with regard to the corrections to wave propagation (through $a_{i}$) and to the interaction between solid and liquid phase (through $\hat{a}_{i}$).

4 Diffusivity estimate

Knowing the crossflow ${\bf U}_{l}-{\bf U}_{s}$ allows to estimate the medium diffusivity in terms of the characteristic scale $L$ of the crystals

[TABLE]

A difficulty arises however because of the oscillatory character of ${\bf U}_{l}-{\bf U}_{s}$. In creeping flow conditions, reversibility of Stokes flows would lead, in the absence of mechanical interactions of the crystals, to pure oscillatory motion of the crystals and of the interstitial liquid; tracer trajectories would be closed and no diffusion would be possible. In order to have diffusion, the effect of inertia and of mechanical interactions must be taken into account, with chaos in the tracer trajectories possibly helping in the process. In this regards, arguments in [21] suggest that tracer excursions in the solid matrix of the order of just a few crystal spacings could be enough to produce trajectory decorrelation. If a condition of large tracer excursion compared to the crystal spacing is satisfied, $\lambda\approx\omega^{-1}|{\bf U}_{l}-{\bf U}_{s}|\gg L$, diffusive behaviours are therefore expected, and Eq. (39) applies.

Considering waves of normalized amplitude $\epsilon$, substitution of Eq. (38) into Eq. (39) gives, after making $L$ dimensionless by rescaling $k_{\omega}L\to L$,

[TABLE]

The large excursion condition $\lambda\gg L$ becomes, using again Eq. (38),

[TABLE]

The drag coefficient in Eqs. (40) and (41) can be estimated as the product $\Gamma\approx nF$ of the drag force on an individual particle $F\approx\mu_{l}L$ and the particle density $n\approx L^{-3}$:

[TABLE]

The condition $\lambda\gg L$, Eq. (41), becomes therefore

[TABLE]

and the diffusivity in Eq. (40), back to dimensional units, takes the form

[TABLE]

The amplification factor with respect to the shear-induced part of the diffusion coefficient in Eq. (2) is ${\cal R}/\epsilon^{2}$. The diffusivity is independent of the volume fraction $\bar{C}$. A more precise calculation carried out in the Appendix, in the case of a dilute fiber suspension, confirms the result within logarithmic corrections.

The estimate in Eq. (44) can be applied to the case of gravity waves in grease ice. In the case of 100 m long waves, 1 mm ice crystals, and an ice viscosity $\nu_{g}\approx\nu_{grease}\approx 0.01\ {\rm m^{2}/s}$, corresponding in dimensionless variables to $\nu_{g}\approx 5\times 10^{-5}$ and $\Gamma\approx 2$, Eqs. (43) and (44) would give

[TABLE]

Using the same ice parameters in the case of 10 m waves, corresponding in dimensionless variables to $\nu_{g}\approx 1.5\times 10^{-3}$ and $\Gamma\approx 0.6$, would give instead ${\cal R}\approx 16$ and $\kappa^{wave}\approx 4\times 10^{-5}\ {\rm m^{2}/s}$.

Note that the diffusivity increment described in Eq. (45) is localized in the viscous boundary layer of thickness $\lambda_{\alpha}=\sqrt{\nu_{g}/\omega}$ near the water surface, which, as already remarked, coincides in most situations of interest with the thickness of the grease ice layer.

5 Concluding remarks

The analysis carried out in the present paper predicts a mechanical contribution to heat and mass diffusion in grease ice stirred by a gravity wave, that exceeds by several orders of magnitude the contribution by shear-induced diffusion. The corresponding increment over the molecular diffusivity is $O(10^{5})$ in the case of salinity, $O(10^{2})$ in the case of heat. Analysis in the Appendix suggests that the contribution of prefactors not included in the analysis may modify the estimates by $O(10^{-2})$ factors, thus it is not possible to conclude that heat diffusion is strongly affected. Less likely that the prefactor contribution is so huge to lead to similar conclusions in the case of salinity. The predicted enhancement of haline diffusion may lead to a shift in the onset of haline convection during ice formation [22].

The mechanism of diffusion described in this paper is not confined to suspension stirring by gravity waves. The mechanism is expected to be active in all situations in which pressure gradients in the flow and non-hydrodynamic interactions of the particles are important.

A question which arises in the case of periodic flows is how much does Stokes flow reversibility in creeping flow conditions affect the development of diffusion. A similar question has attracted recently great attention in the case of shear-induced particle diffusion in periodic shear flows [23]. In that case, irreversibility is restored by particle collisions, provided the shear is above a concentration dependent threshold [24]. In the present case, it is reasonable to expect that a combination of chaos in the tracer trajectories, particle inertia, and dislocations induced by the flow, all play a role in producing the necessary irreversibility for tracer diffusion. The observation is strengthened by the predicted dependence of the diffusivity on the tracer excursion in a wave period in the solid bed compared to the crystal separation. Whether or not some threshold for irreversibility exists, analogous to the one in the case of shear-induced diffusion, remains to be ascertained.

Aknowledgements: This research was supported by FP7 EU project ICE-ARC (Grant agreement No. 603887).

Appendix A Diffusion in a dilute fiber bed

Transverse diffusion in flows in fixed particle beds is possible only if the particles are non-spherical, or if $O(\bar{C}^{2})$ effects are taken into account. The only analytical results available in the first case are those for flows in dilute fiber beds [10]. The case of a random fiber network is considered here, with contact among fibers providing the necessary source of non-hydrodynamic interactions. The length $L$ now indicates the distance along a fiber between different points of contact with other fibers. The volume fraction is expressed in terms of the fiber radius $r$ as $\bar{C}\approx(r/L)^{2}$.

The drag on the liquid phase is obtained, in the dilute limit, from that of an individual fiber, by using slender body theory [25, 26]. In the case of random fiber orientation,

[TABLE]

In order for some form of mechanical diffusion to be present, it is necessary that the displacement $\lambda$ of a tracer relative to the network in a wave period greatly exceeds the correlation length of the velocity fluctuations in the medium. This correlation length can be identified with the Brinkman length [26],

[TABLE]

that is the length above which the velocity perturbation around a fiber is exponentially damped by the image field of the other fibers. In the case of a dilute network, $\lambda_{B}>L$. Redefining ${\cal R}=\lambda/\lambda_{B}$, and exploiting Eqs. (38) and (46) with $\phi_{p}=O(\epsilon)$, gives

[TABLE]

The transverse diffusivity of the medium is given in the large Peclet number by Eq. (31b) in [10]:

[TABLE]

Exploiting Eqs. (38) and (47) gives

[TABLE]

corresponding to a diffusivity reduction with respect to the prediction in Eq. (44), by a factor $\approx 0.02\,|\ln\bar{C}|^{1/2}$.