## 1 Introduction

Jump conditions to be applied across an interface within a viscous fluid have been explored for several decades (Anderson et al., 1998, 2007; Aris, 1962; Delhaye, 1974; Dziubek, 2011; Ishii and Hibiki, 2011; Joseph and Renardy, 1993; Slattery et al., 2007). A general form of the traction jump (the traction is **σn** where **σ** is the stress tensor and **n** a unit vector normal to the interface) involves a flux of momentum across the interface, a possibly anisotropic surface tension and terms including an interface mass density. In pratice, the interface is often supposed to have no mass and the traction is known to undergo a jump especially in two cases: in a shock wave, where the flux of momentum across the interface equals the jump of pressure; and in the presence of surface tension defined as a capillary action due to intermolecular forces at the interface between two immiscible fluids.

Here, we put aside the shock wave and the intrinsic surface tension contributions. In this case, the traction vector is usually supposed to be continuous, for example across phase changes (Hutter and Johnk, 2004). On the contrary, in this paper we show: first, that when a viscous fluid with low Reynolds number crosses an interface with a density jump the traction undergoes a jump; second that this jump takes the mathematical form of an isotropic surface tension that we obtain as a function of the fluid parameters.

The question of a discontinuous traction was first addressed in the context of solid Earth geophysics (Corrieu et al., 1995). Indeed, on a time scale of millions of years, Earth's mantle behaves as a viscous quasi-static fluid that enables convection due to cooling at the external surface and to internal heat sources (decay of radioactive elements). Transported with this movement, Earth's minerals undergo various phase changes due to the large increases in pressure and temperature with depth (and modestly, their variations with latitude and longitude). In Earth, the major phase transitions are located at 410 and 670 km depths and are associated with ∼ 10% density jumps. The physics of these phase changes are complicated as they occur across regions of partial change including a mixture of both phases but each is sharp enough to be modelled as a discontinuity in fluid dynamic simulations of mantle convection. Until now, in simulations the traction jump was supposed to be null at these interfaces (Alboussière et al., 2010; Monnereau et al., 2010; Roberts et al., 2007; Schubert et al., 2001; Sotin and Parmentier, 1989). Concerning the deepest discontinuity in the Earth, a translational convection of the inner solid core within the fluid metallic core has been predicted (Alboussière et al., 2010; Monnereau et al., 2010). This mode relies on a mass flux at the interface, the inner core permanently cristallising on one side and melting on the other. The plausibility of this motion depends however on the traction conditions at the surface of the inner core.

In the context of Earth's convection, Corrieu et al. (1995) noticed however that a stress continuity across a phase jump led to inconsistencies: applying a null traction jump did not yield the same results as solving numerically the Stokes equations within the phase change discontinuity. The authors derived the corresponding jump condition in the special case of a spherical interface, with linear variations of density and viscosity with depth within the transition zone. Ricard (2007) addressed this question again, confirmed that the traction is not continuous and expressed its jump in the case of a plane interface and a uniform viscosity.

The purpose of the present paper is to examine this question in a more general framework, for interfaces of arbitrary curvature, density and Newtonian viscosity variations. The model is justified in the geophysical context but can be considered as rather specific in general fluid mechanics: it is fully dominated by viscous stresses (small Reynolds number) and the density variations of the background are fixed and represented by a discontinuity. Notice however that the problem appearing in the case of a creeping flow should also be present in the case of finite Reynolds numbers. We show that in the case of a density jump there must be an equivalent surface tension to ensure global conservation of momentum, and that this tension is given in terms of the mass transfer flux and of the density and viscosity changes by Eqs. (25)–(26). Upon applications, this effect may influence the inner core crystallisation.

In the next sections, we point out how the traction jump appears in two simple cases. Then, the following section is devoted to the expression and proof of the traction jump in a general case: we look for weak solutions of the Stokes equations that are smooth on either side of an arbitrarily curved interface, and thus obtain the announced jump conditions.

## 2 Existence and necessity of a traction jump

Searching the interface condition is related to the following question: if a velocity and stress field is a solution of the Stokes equations for a continuous but rapidly varying density, can this solution be approximated by the one for a sharp interface and what are the corresponding jump conditions? We illustrate the problem with two simple situations, first with a flat interface, then with a curved interface.

### 2.1 Existence of a traction jump across a flat interface

We consider a 2D compressible, viscous, quasi-static liquid, in steady-state regime, through the rectangular volume −h < z < h and −π/k < x < π/k. A mass flux ψ_{0} cos kx is injected at z = −h and extracted at z = h. To mimic a phase transition that would occur near z = 0, we assume that the density ρ(z) changes continuously from ρ^{−} to ρ^{+} between z = −ɛ and z = ɛ and we solve for the flow **υ** (x, y) in the whole domain. Considering that a phase change occurs across a very thin z-interval, this situation is expected to close to the situation where the density changes discontinuously at z = 0 and where jump conditions are implemented across the discontinuity. We show however that the continuous approach when ɛ→0 does not converge to the discontinuous approach where the traction is supposed continuous.

The Stokes equations governing a viscous linear quasi-static fluid are:

$$\nabla \cdot \left(\rho \upsilon \right)=0\text{,}$$ | (1) |

$$\nabla \cdot \sigma =0\text{,}$$ | (2) |

$$\sigma =-pI\hspace{0.17em}+\eta \left(\nabla \upsilon +{\left(\nabla \upsilon \right)}^{T}\right)+\lambda \left(\nabla \cdot \upsilon \right)I\text{,}$$ | (3) |

**υ**its velocity,

**σ**the stress tensor, p the pressure, η and λ the shear and bulk viscosities, and

**I**the identity tensor.

According to the first relation, ρ**υ** can be sought of the form $\rho {\upsilon}_{x}=-{\partial}_{z}\Psi ,\text{\hspace{0.28em}}\rho {\upsilon}_{z}=+{\partial}_{x}\Psi $; with the curl of relation (2), and assuming constant η, it implies that

$$\Delta \left(\nabla \cdot \left(\frac{1}{\rho}\nabla \mathrm{\Psi}\right)\right)=0\text{.}$$ | (4) |

At z = ± h we impose the sinusoidal vertical flux:

$${\partial}_{z}\mathrm{\Psi}\left(x,z=\pm h\right)=0\text{,}$$ | (5) |

$${\partial}_{x}\mathrm{\Psi}\left(x,z=\pm h\right)={\Psi}_{0}\text{\hspace{0.28em}}\mathrm{cos}\text{\hspace{0.17em}}kx\text{.}$$ | (6) |

Solutions of the form $\Psi =\stackrel{\u02c6}{\mathrm{\Psi}}$ (z) sin (kx) are easy to find numerically as $\stackrel{\u02c6}{\Psi}$ is the solution of a fourth order linear and homogeneous differential system in z. Using a Runge–Kutta method, four elementary solutions are computed in each domain where the parameters are continuous (i.e., in the whole domain when the density evolves continuously, in the two subdomains z < 0 and z > 0 when the density is discontinuous), then the general solution is found by linearly combining these solutions, with coefficients determined by solving the linear system corresponding to the limit and jump conditions. In that way we compute three types of solutions:

- • a solution with a smooth but rapidly varying density ρ = ρ
_{ɛ}(Fig. 1, solid black line). We choose ρ_{ɛ}(z) = ρ^{–}+ (ρ^{+}– ρ^{–}) (tanh (z/ɛ)+ 1)/2, which, for small ɛ is close to the function ${\rho}_{0}={\rho}^{-}{1}^{-}+{\rho}^{+}{1}^{+}\text{\hspace{0.28em}}\text{where}\text{\hspace{0.28em}}{1}^{\pm}$ denotes the characteristic functions of the half-spaces $z\gtrless 0$. - • a solution with a discontinuous density ρ = ρ
_{0}(Fig. 1, dashed red line), a continuous mass flux (ρυ_{z}) and traction at z = 0. These are the usual conditions for a flow crossing a phase change, neglecting surface tension. - • a solution with a discontinuous density ρ = ρ
_{0}, a continuous mass flux and a discontinuous traction at z = 0 as proposed later in this paper.

Although the vertical velocities are visually similar (Fig. 2, top row), the results for the horizontal velocity (bottom row) are surprisingly different. This suggests that the limit of a diffuse layer (solution [a], left panel) is not a sharp interface with a continuous traction (solution [b], middle panel) but a sharp interface with an equivalent dynamic surface tension (solution [c], right panel).

To understand where the problem lies, we plot on Fig. 3, the profiles of the shear stress σ_{xz}(x = π/(2k), z) (left) and vertical stress σ_{zz}(x = 0, z) (right), at x = 0. It is obvious that the shear stress profile computed for the rapid and continuous density change (solution [a], solid line, left panel) looks rather discontinuous and can hardly converge toward the smooth solution computed with the usual jump conditions on interfaces (solution [b], blue dotted-dashed line, left panel). The profiles for the vertical stress are continuous but remain clearly different (right panel). The solution computed with a continuous mass flux and a discontinuous traction at z = 0 as proposed later in this paper (solution [c] red dashed lines) provides clearly a much better fit to the stresses computed with a continuous density change (solution [a], solid lines).

This simple example, which can be interpreted as a naive view of mantle convection with a hot mantle crossing a phase transition in between two subducting slabs, illustrates therefore a situation where the usual jump conditions on interface do not hold and where the shear stress is discontinuous on an interface without intrinsic surface tension.

### 2.2 The necessity of a traction jump

Let us now explain in a heuristic way why the traction must have a jump across a density jump with a mass transfer. As in the example before, the interface is considered as a layer in between z = −ɛ and z = +ɛ inside which Eqs. (1)–(3) hold and all the fluid parameters vary significantly faster in the z-direction than in the x-direction.

The thought experiment consisting in considering ɛ→0 is called in most textbooks the “pillbox” argument. Mathematically speaking, we start from a family of rapidly varying but still smooth densities ρ_{ɛ}, and viscosities λ_{ɛ} and η_{ɛ}, converging pointwise when ɛ goes to zero to some functions ρ, λ and η that are possibly discontinuous at z = 0 – which means that the limiting sharp interface is located at z = 0. Then we assume that with those “mollified” coefficients ρ_{ɛ}, λ_{ɛ}, η_{ɛ}, the 2D Stokes equations have solutions **υ**_{ɛ} converging pointwise to **υ**, and we look for jump conditions on **υ**. The general idea is that ρ_{ɛ}, λ_{ɛ}, η_{ɛ}, **υ**_{ɛ}, and their x-derivatives remain bounded in the limit ɛ→0, whereas the z-derivative of any of those quantities that is discontinuous in the limit blows up as 1/ɛ. For simplicity, we omit the subscript ɛ in what follows.

In 2D, the Stokes Eqs. (1)–(3) read:

$${\partial}_{z}\left(\rho {\upsilon}_{z}\right)+{\partial}_{x}\left(\rho {\upsilon}_{x}\right)=0$$ | (7) |

$${\partial}_{z}{\sigma}_{zz}+{\partial}_{x}{\sigma}_{xz}=0$$ | (8) |

$${\partial}_{z}{\sigma}_{xz}+{\partial}_{x}{\sigma}_{xx}=0$$ | (9) |

$${\sigma}_{xz}=\eta \left({\partial}_{x}{\upsilon}_{z}+{\partial}_{z}{\upsilon}_{x}\right)$$ | (10) |

$${\sigma}_{zz}=-p+2\eta {\partial}_{z}{\upsilon}_{z}+\lambda \left({\partial}_{x}{\upsilon}_{x}+{\partial}_{z}{\upsilon}_{z}\right)$$ | (11) |

$${\sigma}_{xx}=-p+2\eta {\partial}_{x}{\upsilon}_{x}+\lambda \left({\partial}_{x}{\upsilon}_{x}+{\partial}_{z}{\upsilon}_{z}\right)\text{.}$$ | (12) |

Eq. (7) readily shows that ρυ_{z} must remain continuous in the limit ɛ→0 (otherwise, ∂_{x}(ρυ_{x}) = −∂_{z}(ρυ_{z}) would blow up). If we denote by $\u301a\text{\hspace{0.28em}}\u301b$ the jump of a quantity across the interface, this usual jump condition reads $\u301a\rho {\upsilon}_{z}\u301b=0$. By the same argument, the usual jump conditions are obtained by considering that equations (8) and (9) imply $\u301a{\sigma}_{zz}\u301b=0\text{\hspace{0.28em}}\text{and}\text{\hspace{0.28em}}\u301a{\sigma}_{xx}\u301b=0$, respectively, provided that ∂_{x}σ_{xz} and ∂_{x}σ_{xx} be bounded. However, the latter assumption is unlikely to hold true if the mass transfer flux ρυ_{z} is nonzero. Indeed, subtracting (11) to (12) we obtain

$${\sigma}_{xx}={\sigma}_{zz}+2\eta \left({\partial}_{x}{\upsilon}_{x}-{\partial}_{z}{\upsilon}_{z}\right)\text{.}$$ | (13) |

The important point is that the vertical velocity becomes discontinuous when the thickness ɛ→0 and that ∂_{z}υ_{z} blows up. As suggested by Fig. 3, σ_{zz} remains bounded (and continuous) and as ∂_{x}υ_{x} is also bounded, σ_{xx} must blow up too. In (9), the two terms blow up as 1/ɛ and one cannot infer that $\u301a{\sigma}_{xz}\u301b=0$ anymore.

To derive the appropriate jump, we use (9) and (13) and find that

$${\partial}_{z}{\sigma}_{xz}=-{\partial}_{x}{\sigma}_{xx}=-{\partial}_{x}{\sigma}_{zz}-2{\partial}_{x}\left(\eta \left({\partial}_{x}{\upsilon}_{x}-{\partial}_{z}{\upsilon}_{z}\right)\right)$$ | (14) |

In the right-hand side, only the term in ∂_{z}υ_{z} is singular in the limit ɛ→0, thus by integration across the layer (Anderson et al., 1998) on a height 2δ such that $\epsilon \ll \delta \ll 1$ we get $\u301a{\sigma}_{xz}\u301b=2{\mathrm{lim}}_{\epsilon \to 0}{\partial}_{x}{\int}_{-\delta}^{\delta}\eta {\partial}_{z}{\upsilon}_{z}\text{\hspace{0.28em}}dz\text{.}$ With the mass conservation (7) it implies $\u301a{\sigma}_{xz}\u301b=-2{\mathrm{lim}}_{\epsilon \to 0}{\partial}_{x}{\int}_{-\delta}^{\delta}\eta \left({\partial}_{z}\rho \right){\upsilon}_{z}/\rho \text{\hspace{0.28em}}\text{d}z$ or since ρυ_{z} is continuous:

$$\u301a{\sigma}_{xz}\u301b=-2{\partial}_{x}\left(\rho {v}_{z}\underset{\epsilon \to 0}{{\displaystyle \mathrm{lim}}}{\int}_{-\delta}^{\delta}\eta \left({\partial}_{z}\rho \right)/{\rho}^{2}\text{d}z\right)\text{.}$$ | (15) |

This can be written

$$\u301a{\sigma}_{xz}\u301b=2{\partial}_{x}\left(\u301a\nu \u301b\rho {\upsilon}_{z}\right)$$ | (16) |

$$\u301a\nu \u301b=\underset{\epsilon \to 0}{{\displaystyle \mathrm{lim}}}{\int}_{-\delta}^{\delta}\eta {\partial}_{z}(1/\rho )\text{d}z$$ | (17) |

The jump conditions for a plane interface are therefore:

$$\begin{array}{l}\u301a\rho {\upsilon}_{z}\u301b=0\text{\hspace{1em}}\text{\hspace{1em}}\u301a{\sigma}_{zz}\u301b=0\\ \u301a{\upsilon}_{x}\u301b=0\text{\hspace{1em}}\text{\hspace{1em}}\u301a{\sigma}_{xz}\u301b=2{\partial}_{x}\left(\u301av\u301b\rho {\upsilon}_{Z}\right).\end{array}$$ | (18) |

With a constant viscosity within the interface we get $\u301a\nu \u301b=\eta \u301a1/\rho \u301b$, and since the first relation implies $\u301a1/\rho \u301b\rho {\upsilon}_{z}=\u301a{\upsilon}_{z}\u301b$ the fourth equation also reads $\u301a{\sigma}_{xz}\u301b=2{\partial}_{x}\left(\eta \u301a{\upsilon}_{z}\u301b\right)\text{.}$ The correct jump condition yields $\u301a{\partial}_{z}{\upsilon}_{x}\u301b-{\partial}_{x}\u301a{\upsilon}_{z}\u301b=0$ while the traditional condition $\u301a{\sigma}_{xz}\u301b=0$ implies $\u301a{\partial}_{z}{\upsilon}_{x}\u301b+{\partial}_{x}\u301a{\upsilon}_{z}\u301b=0!$

According to (13), σ_{xx} becomes infinite in the limit ɛ→0. Using the same kind of integration as for the proof of (16), its integral across the infinitely thin interface is

$$s\u2254\underset{\epsilon \to 0}{{\displaystyle \mathrm{lim}}}{\int}_{-\delta}^{\delta}{\sigma}_{xx}\text{d}z=-2\u301a\nu \u301b\rho {\upsilon}_{z}\text{.}$$ | (19) |

When the flow crosses a density interface, the pillbox argument cannot be used without including lateral terms because, in the same way as in the presence of intrinsic surface tension (Anderson et al., 1998), the lateral stress diverges when the pillbox shrinks and thus the lateral integrals do not vanish in the zero-height limit.

To summarize, we expect that for a viscous fluid (η ≠ 0) with mass transfer at the interface (ρυ_{z} ≠ 0), both p and σ_{xx} become singular – like across an interface endowed with intrinsic surface tension – while ρυ_{z} and υ_{x} are continuous, and the other stress components σ_{xz} and σ_{zz} may have jumps (for a flat interface, the latter is not affected by surface tension).

As we have seen that a rapid density change implies the existence of an equivalent surface tension, we should now verify that the normal stress discontinuity occurring across a curved interface (the Young–Laplace condition) is also valid in the case of a rapid density change.

### 2.3 Existence of a traction jump across a curved interface

We consider a simple flow diverging from a point source with radial velocity $\upsilon =\upsilon \left(r\right)\stackrel{\u02c6}{r}$, with a density ρ (r) varying rapidly from ρ^{−} to ρ^{+} across a shell located between the radii R − δ and R + δ. In this spherically symmetric situation the solution of Stokes equations is even more surprising since we can prove analytically (and very simply) the existence of a dynamic surface tension. Indeed, with η = η (r) and λ = λ (r), the only nonzero stress is ${\sigma}_{rr}=-p+2\eta {\partial}_{r}\upsilon +\lambda {\partial}_{r}\left({r}^{2}\upsilon \right)/{r}^{2}$ and ${\sigma}_{\theta \theta}={\sigma}_{\varphi \varphi}=-p+2\eta \upsilon /r+\lambda {\partial}_{r}\left({r}^{2}\upsilon \right)/{r}^{2}$ and the equilibrium equation yields

$${\partial}_{r}{\sigma}_{rr}+4\eta {\partial}_{r}\left(\frac{\upsilon}{r}\right)=0\text{.}$$ | (20) |

The mass conservation equation implies

$$\rho \upsilon =\frac{C}{{r}^{2}}$$ | (21) |

$${\partial}_{r}{\sigma}_{rr}+4\eta \frac{\rho \upsilon}{r}{\partial}_{r}\left(1/\rho \right)-12\eta \frac{\upsilon}{r}=0\text{.}$$ | (22) |

The last term is bounded, ρυ/r is continuous, therefore by integration across the discontinuity

$$\u301a{\sigma}_{rr}\u301b=-4\frac{\rho \upsilon}{R}\underset{\epsilon \to 0}{{\displaystyle \mathrm{lim}}}{\int}_{R-\delta}^{R+\delta}\eta {\partial}_{r}\left(1/\rho \right)\text{d}r=\frac{2s}{R}\text{.}$$ | (23) |

As in the Young–Laplace law, the jump of normal stress is not zero but proportional to the total curvature 2/R and to the surface tension which is here $-2\u301a\nu \u301b\rho \upsilon $, in agreement with what we found across a flat interface. This effect may be large, as the jump $\left(2s/R\sim 4\eta \u301a1/\rho \u301b\rho \upsilon /R\right)$ and the viscous stress $\left(\sim \eta \upsilon /R\right)$ have similar magnitudes. In this example, that can be seen as a naive view of the crystallising inner core, the usual jump conditions on interface also do not hold, and the normal stress is discontinuous on an interface without intrinsic surface tension.

## 3 Jump conditions

### 3.1 Statement of jump conditions

We now give the general jump conditions for a general 3D fluid having a curved interface. We consider a fluid verifying eqs. (1)–(3) everywhere, and undergoing a density jump when crossing a static surface $\Sigma =\left\{x\left|{q}^{3}=0\right.\right\}$, with q^{3} a local coordinate normal to the surface. We suppose that this mathematical surface corresponds physically to a very small layer in which the Stokes equations remain valid. We neglect intrinsic surface tension on the interface Σ. The jump of any piecewise-smooth quantity f across the interface at $x\in \Sigma $ is defined as $\u301af\u301b={f}^{+}-{f}^{-}\u2254{\mathrm{lim}}_{\epsilon \to {0}^{+}}f\left(x+\epsilon n\right)-{\mathrm{lim}}_{\epsilon \to {0}^{+}}f\left(x-\epsilon n\right)$ where **n** denotes the unit normal to the interface at point **x** and pointing towards positive q^{3}. In addition, we use the subscript _{T} for the tangential part ${\upsilon}_{T}=\upsilon -\left(\upsilon \cdot n\right)n$ of a vector field υ on Σ. Accordingly, ${\nabla}_{T}$ and ${\nabla}_{T}\cdot $ stand for the tangential gradient and tangential divergence on Σ (Slattery et al., 2007). Then in particular, $\kappa ={\nabla}_{T}\cdot n$ is the sum of principal curvatures of Σ. The main result of the paper is that, apart from the usual jump conditions

$$\u301a\rho \upsilon \cdot n\u301b=0,\text{\hspace{1em}}\u301a{\upsilon}_{T}\u301b=0\text{,}$$ | (24) |

$$\u301a\sigma n\u301b=s\kappa n-{\nabla}_{T}s\text{,}$$ | (25) |

$$s\u2254-2\u301a\nu \u301b\rho \upsilon \cdot n$$ | (26) |

$$\u301a\frac{1}{\rho}\u301b\text{\hspace{0.28em}}\rho \upsilon \cdot n=\u301a\upsilon \cdot n\u301b$$ | (27) |

The usual jump conditions on a phase change interface consist of eq. (24) together with $\u301a\sigma n\u301b=0$ (Hutter and Johnk, 2004). The latter is recovered from (25) when the dynamic surface tension s vanishes identically, which is the case for an impermeable interface (thus with null mass flux) or if the density is uniform. The traction continuity is recovered too when the interface is flat and the dynamic surface tension s is uniform.

For the sake of simplicity, inertia, surface forces, surface mass density and any a priori intrinsic surface tension have been neglected. These effects are taken into account by well-known jump conditions (Slattery et al., 2007). The main result of the paper is that without any intrinsic surface tension there still exists another surface tension, but of a dynamic origin, which is explicitely given as a function of the fluid parameters by (26).

### 3.2 Proof of jump conditions

Finding the jump conditions amounts to investigate the (very) weak solutions to the Stokes equations which are expected from a passage to the limit ɛ→0 as in Section 2.2. For that purpose we look for jump conditions on solutions of (1)–(3) in the sense of distributions, with discontinuous coefficients ρ, λ, η across a smooth surface Σ.

We denote by ψ(**x**) the signed distance of the point **x** to the surface. The equation of the surface is ψ(**x**) = 0, and in a vicinity of the surface ψ satisfies the eikonal equation $\u2225\nabla \psi \left(x\right)\u2225=1$ and defines the unit normal $n\u2254\nabla \psi \left(x\right)$. We denote by x_{1}, x_{2}, x_{3} = x, y, z the Cartesian coordinates of **x**, and we adopt the usual convention of summation over repeated indices. For convenience we rewrite the eq. (2) with indices:

$${\partial}_{j}\left(-p+\lambda {\partial}_{i}{\upsilon}_{i}\right)+{\partial}_{i}\left(\eta \left({\partial}_{i}{\upsilon}_{j}+{\partial}_{j}{\upsilon}_{i}\right)\right)=0,\text{\hspace{1em}}j\in \left\{\mathrm{1,2,3}\right\}\text{,}$$ | (28) |

**L**(p, υ) = 0 where

**L**is the Stokes operator

$$\left(p,\upsilon \right)\mapsto L\left(p,\upsilon \right)\u2254\nabla \left(-p+\lambda \nabla \cdot \upsilon \right)+\nabla \cdot \left(\eta \left(\nabla \upsilon +{\left(\nabla \upsilon \right)}^{\text{T}}\right)\right)\text{.}$$ | (29) |

A usual method to get jump conditions is to find the discontinuous solutions of the corresponding operator (Bedeaux et al., 1976). Thus, we assume that the velocity field is of the form

$$\upsilon ={\upsilon}^{-}{1}^{-}+{\upsilon}^{+}{1}^{+}\text{,}$$ | (30) |

$$p={p}^{-}{1}^{-}+{p}^{+}{1}^{+}+P\delta \text{,}$$ | (31) |

**υ**

^{±}, p

^{±}are smooth functions of

**x**on the half-spaces

$$\left\{x=\left({x}_{1},{x}_{2},{x}_{3}\right);\psi \left(x\right)\ge 0\right\},\left\{x=\left({x}_{1},{x}_{2},{x}_{3}\right);\psi \left(x\right)\le 0\right\}\text{,}$$ | (32) |

**x**on Σ and its vicinity.

In what follows, we repeatedly use the following calculus formulas (Bedeaux et al., 1976; Gel’fand and Shilov, 1968, chap. 3.1):

$${\partial}_{j}1\pm =\pm \left({\partial}_{j}\psi \right)\delta ,\text{\hspace{1em}}{\partial}_{j}\delta =\left({\partial}_{j}\psi \right){\delta}^{\prime}\text{,}$$ | (33) |

Then, the derivatives of p and **υ** in the sense of distributions are

$${\partial}_{j}p=\left({\partial}_{j}{p}^{-}\right){1}^{-}+\left({\partial}_{j}{p}^{+}\right){1}^{+}+\u301ap\u301b\left({\partial}_{j}\psi \right)\delta +P\left({\partial}_{j}\psi \right){\delta}^{\prime}\text{,}$$ | (34) |

$${\partial}_{j}\upsilon =\left({\partial}_{j}{\upsilon}^{-}\right){1}^{-}+\left({\partial}_{j}{\upsilon}^{+}\right){1}^{+}+\u301a\upsilon \u301b\left({\partial}_{j}\psi \right)\delta \text{.}$$ | (35) |

For convenience, we extend **υ**^{±} in a smooth manner to the other sides of Σ, so that the jump $\u301a\upsilon \u301b$ is extended to a smooth function in the vicinity of Σ, merely defined as $\u301a\upsilon \u301b={\upsilon}^{+}-{\upsilon}^{-}$. In the same way we extend p^{±} and $P$.

For the sake of clarity, we start with λ and η continuous across the interface and postpone our more elaborate calculation intended for general interfaces to the next section.

#### 3.2.1 Continuous rheology

Let us first derive jump conditions in the special case when all products in (28) are well defined, that is, assuming that both viscosities λ and η are regular, and more especially constant across Σ.

Differentiating once more in the sense of distributions and using ${\partial}_{j}\psi ={n}_{j}$ and (33) we find that the Stokes equations **L**(p, **υ**) = 0 are equivalent to **L**(p^{±}, **υ**^{±}) = 0 and **m**(p, **υ**) = 0, where **m**(p, **υ**) is the vector-valued distribution supported by Σ whose components are

$${m}_{j}\left(p,\mathrm{\upsilon}\right)=\begin{array}{lll}-\u301ap\u301b{n}_{j}\delta \hfill & -\left({\partial}_{j}P\right)\delta \hfill & -P{n}_{j}{\delta}^{\prime}\hfill \\ +\u301a\lambda {\partial}_{i}{\upsilon}_{i}\u301b{n}_{j}\delta \hfill & +{\partial}_{j}\left(\u301a\lambda {\upsilon}_{i}\u301b{n}_{i}\right)\delta \hfill & +\u301a\lambda {\upsilon}_{i}\u301b{n}_{i}{n}_{j}{\delta}^{\prime}\hfill \\ +\u301a\eta {\partial}_{i}{\upsilon}_{j}\u301b{n}_{i}\delta \hfill & +{\partial}_{i}\left(\u301a\eta {\upsilon}_{j}\u301b{n}_{i}\right)\delta \hfill & +\u301a\eta {\upsilon}_{j}\u301b{n}_{i}{n}_{i}{\delta}^{\prime}\hfill \\ +\u301a\eta {\partial}_{j}{\upsilon}_{i}\u301b{n}_{i}\delta \hfill & +{\partial}_{i}\left(\u301a\eta {\upsilon}_{i}\u301b{n}_{j}\right)\delta \hfill & +\u301a\eta {\upsilon}_{i}\u301b{n}_{j}{n}_{i}{\delta}^{\prime},\hfill \end{array}$$ | (36) |

Now we use that for any smooth functions a(**x**) and b(**x**), the equality aδ + bδ′ =0 is equivalent to b = 0 and $a={n}_{j}{\partial}_{j}b$ on the interface Σ (here, ${\delta}^{\prime}={n}_{j}{\partial}_{j}\delta $ is a normal derivative, eq. (33), so that $a\delta +b{\delta}^{\prime}$ applied to any test function f, is the value of $af-{\partial}_{j}\left({n}_{j}bf\right)=\left(a-b{\partial}_{j}{n}_{j}-{n}_{j}{\partial}_{j}b\right)f-b{n}_{j}{\partial}_{j}f$ taken on the interface). In this paper, b will represent functions defined on the interface only, thus their extension outside Σ does not play any role and we can assume, for convenience and without loss of generality, that their normal derivative is zero. Thus, $a\delta +b{\delta}^{\prime}=0$ will be equivalent to b = 0 and a = 0 on Σ (Bedeaux et al., 1976, appendix A).

Thus, collecting the δ′ terms of (36) and equating their sum to zero, we infer that

$$P{n}_{j}=\u301a\lambda \upsilon \cdot n\u301b\text{\hspace{0.28em}}{n}_{j}+\u301a\eta {\upsilon}_{j}\u301b+\u301a\eta \upsilon \cdot n\u301b\text{\hspace{0.28em}}{n}_{j}$$ | (37) |

**n**. In other words, the tangential velocity must be continuous when viscosities are so, and then we merely have

$$P=\u301a\left(\lambda +2\eta \right)\upsilon \cdot n\u301b$$ | (38) |

_{j}(p,

**υ**) must cancel out, which leads to

$$\begin{array}{l}{\partial}_{j}P=-\u301ap\u301b{n}_{j}+\u301a\lambda \nabla \cdot \upsilon \u301b{n}_{j}+\u301a\eta ({\partial}_{i}{\upsilon}_{j}+{\partial}_{j}{\upsilon}_{i})\u301b{n}_{i}\\ +{\partial}_{j}\left(\u301a\lambda \upsilon \cdot n\u301b\right)+{\partial}_{i}(\u301a\eta {\upsilon}_{j}\u301b{n}_{i})+{\partial}_{i}\left(\u301a\eta {\upsilon}_{i}\u301b{n}_{j}\right)\end{array}$$ | (39) |

$$={\u301a\sigma n\u301b}_{j}+{\partial}_{j}\left(\u301a\lambda \upsilon \cdot n\u301b\right)+{\partial}_{i}\left(\u301a\eta {\upsilon}_{j}\u301b{n}_{i}\right)+{\partial}_{i}\left(\u301a\eta {\upsilon}_{i}\u301b{n}_{j}\right)\text{.}$$ | (40) |

Eliminating ${\partial}_{j}P$ using (38), we obtain

$${\u301a\sigma n\u301b}_{j}={\partial}_{j}\left(2\u301a\eta \upsilon \cdot n\u301b\right)-{\partial}_{i}\left(\u301a\eta {\upsilon}_{j}\u301b{n}_{i}\right)-{\partial}_{i}\left(\u301a\eta {\upsilon}_{i}\u301b{n}_{j}\right)\text{.}$$ | (41) |

We take the inner product of this equation by the normal **n**. Various simplifications occur as **n** is a unit vector, as the surfaces ψ(**x**) = cst are parallel and as $\u301a\upsilon \u301b$ is parallel to **n**:

$${n}_{i}{n}_{i}=1,\text{\hspace{0.28em}}{n}_{i}{\partial}_{i}{n}_{j}={n}_{i}{\partial}_{i}{\partial}_{j}\psi ={n}_{i}{\partial}_{j}{n}_{i}={\partial}_{j}\left({n}_{i}{n}_{i}\right)/2=0,\text{\hspace{0.28em}}\u301a{\upsilon}_{i}\u301b=\u301a\upsilon \cdot n\u301b{n}_{i}\text{,}$$ | (42) |

$$\u301a\left(\sigma n\right)\cdot n\u301b=-2\u301a\eta \upsilon \cdot n\u301b\kappa $$ | (43) |

$$\kappa \u2254{\partial}_{i}{n}_{i}=\Delta \psi \text{.}$$ | (44) |

Now, taking the inner product of (41) by any tangent vector $\text{T}\perp \text{n}$, and using again $\u301a{\upsilon}_{i}\u301b=\u301a\mathrm{\upsilon}\cdot n\u301b{n}_{i}$ we arrive at

$$\u301a\left(\sigma n\right)\cdot T\u301b={\partial}_{j}\left(2\u301a\eta \upsilon \cdot n\u301b\right){T}_{j}-2{\partial}_{i}\left(\u301a\eta \upsilon \cdot n\u301b{n}_{i}{n}_{j}\right){T}_{j}\text{,}$$ | (45) |

_{i}∂

_{i}n

_{j}= 0 and T

_{j}n

_{j}= 0. Hence the other jump condition

$$\u301a\left(\sigma n\right)\cdot T\u301b=2{\partial}_{T}\u301a\eta \upsilon \cdot n\u301b\text{.}$$ | (46) |

Since **T** is any tangent vector, the relations (43)–(46) coincide with (25)–(26) when ν = η/ρ (which is indeed the case when η = cst). Of course (43)–(46) also reduces to (18) for a flat interface (κ = 0) and to (23) for a spherical symmetry (κ = 2/R).

#### 3.2.2 Discontinuous rheology

Let us now extend this computation to possibly discontinuous viscosities. As a consequence, the distributions λ∂_{i}υ_{j} and η∂_{i}υ_{j} in (3) are not well defined (product of a Heaviside function by a Dirac mass at its point of discontinuity).

To overcome this difficulty, we first observe that by the mass eq. (1), we have $\u301a\rho \upsilon \cdot n\u301b=0$ across Σ. This is the most classical, Rankine–Hugoniot jump condition, which can easily be recovered by the method used above for instance. If ρ is genuinely discontinuous, this implies that $\upsilon \cdot n$ also has a discontinuity across Σ, and as a consequence the distribution ${\partial}_{n}\left(\upsilon \cdot n\right)\u2254n\cdot \nabla \left(\upsilon \cdot n\right)={n}_{i}{\partial}_{i}\left({n}_{j}{\upsilon}_{j}\right)={n}_{i}{n}_{j}{\partial}_{i}{\upsilon}_{j}$ involves a Dirac mass. However, we also notice that

$${\partial}_{n}\left(\upsilon \cdot n\right)=\left(1/\rho \right){\partial}_{n}\left(\rho \upsilon \cdot n\right)+\rho \upsilon \cdot n{\partial}_{n}\left(1/\rho \right)\text{,}$$ | (47) |

$$\begin{array}{ccc}\hfill \lambda {\partial}_{n}\left(\upsilon \cdot n\right)\hfill & \hfill =\hfill & \hfill \left(\lambda /\rho \right){\partial}_{n}\left(\rho \upsilon \cdot n\right)+\rho \upsilon \cdot n{\partial}_{n}\mu ,\hfill \\ \hfill \eta {\partial}_{n}\left(\upsilon \cdot n\right)\hfill & \hfill =\hfill & \hfill \left(\eta /\rho \right){\partial}_{n}\left(\rho \upsilon \cdot n\right)+\rho \upsilon \cdot n{\partial}_{n}\nu ,\hfill \end{array}$$ | (48) |

$$\begin{array}{l}\lambda {\partial}_{n}\left(\upsilon \cdot n\right)={\lambda}^{-}{\partial}_{n}\left({\upsilon}^{-}\cdot n\right){1}^{-}+{\lambda}^{+}{\partial}_{n}\left({\upsilon}^{+}\cdot n\right){1}^{+}+\u301a\mu \rho \upsilon \cdot n\u301b\delta ,\\ \eta {\partial}_{n}\left(\upsilon \cdot n\right)={\eta}^{-}{\partial}_{n}\left({\upsilon}^{-}\cdot n\right){1}^{-}+{\eta}^{+}{\partial}_{n}\left({\upsilon}^{+}\cdot n\right){1}^{+}+\u301a\nu p\upsilon \cdot n\u301b\delta ,\end{array}$$ | (49) |

Still, the mass equation does not give any information on the tangential velocity **υ**_{T} so that the corresponding products $\eta {\partial}_{n}{\left({\upsilon}_{T}\right)}_{i}$ in the stress eq. (3) are not well defined either. As we have seen before, the continuity of tangential velocity across Σ is a consequence of the stress equation when η itself is continuous. Otherwise, we have to assume the continuity of tangential velocity in order to deal with the related products in (2). Thus, we assume that the tangential velocity is continuous, that is $\u301a{\upsilon}_{T}\u301b=0$ across Σ to avoid undefined products. In this way, the products λ∂_{j}υ_{i} and η∂_{j}υ_{i} are not singular, except for the component on n_{j}n_{i} that is given by (49). Thus, we deduce that

$$\begin{array}{l}\lambda {\partial}_{j}{\upsilon}_{i}={\lambda}^{-}\left({\partial}_{j}{\upsilon}_{j}^{-}\right){1}^{-}+{\lambda}^{+}\left({\partial}_{j}{\upsilon}_{i}^{+}\right){1}^{+}\\ +\text{\hspace{0.28em}}\u301a\mu \text{\hspace{0.17em}}\rho \upsilon \cdot n\u301b{n}_{j}{n}_{i}\delta ,\\ \eta {\partial}_{j}{\upsilon}_{i}={\eta}^{-}\left({\partial}_{j}{\upsilon}_{j}^{-}\right){1}^{-}+{\eta}^{+}\left({\partial}_{j}{\upsilon}_{i}^{+}\right){1}^{+}\\ +\text{\hspace{0.28em}}\u301a\nu \text{\hspace{0.17em}}\rho \upsilon \cdot n\u301b{n}_{j}{n}_{i}\delta ,\end{array}$$ | (50) |

In case λ = cst, we do not need to introduce μ, and the jump formula above and below are valid with λ/ρ instead of μ [it suffices to apply (35) to $\left(\upsilon \cdot n\right)$ and multiply by λ]. In the same way, if η = cst we can just replace **υ** by η/ρ.

With these relations we find that the vector-valued distribution **m**(p, **υ**) has components

$${m}_{j}\left(p,\upsilon \right)=\begin{array}{lll}-\u301ap\u301b{n}_{j}\delta \hfill & -\left({\partial}_{j}P\right)\delta \hfill & -P{n}_{j}{\delta}^{\prime}\hfill \\ +\u301a\lambda {\partial}_{i}{\upsilon}_{i}\u301b{n}_{j}\delta \hfill & +{\partial}_{j}(\u301a\mu \rho \upsilon \cdot n\u301b)\delta \hfill & +\u301a\mu \rho \upsilon \cdot n\u301b{n}_{j}{\delta}^{\prime}\hfill \\ +\u301a\eta {\partial}_{i}{\upsilon}_{j}\u301b{n}_{i}\delta \hfill & +{\partial}_{i}\left(\u301a\nu \rho \upsilon \cdot n\u301b{n}_{j}{n}_{i}\right)\delta \hfill & +\u301a\nu \rho \upsilon \cdot n\u301b{n}_{j}{\delta}^{\prime}\hfill \\ +\u301a\eta {\partial}_{j}{\upsilon}_{i}\u301b{n}_{i}\delta \hfill & +{\partial}_{i}\left(\u301a\nu \rho \upsilon \cdot n\u301b{n}_{i}{n}_{j}\right)\delta \hfill & +\u301a\nu \rho \upsilon \cdot n\u301b{n}_{j}{\delta}^{\prime}.\hfill \end{array}$$ | (51) |

This expression corresponds to the continuous case if we replace μ and **υ** with λ/ρ and η/ρ and use (42). Now, **m**(p, **υ**) = 0 readily implies (collecting δ′) that

$$P=\u301a\left(\mu +2\nu \right)\rho \upsilon \cdot n\u301b\text{,}$$ | (52) |

$${\partial}_{j}P={\u301a\sigma n\u301b}_{j}+{\partial}_{j}\u301a\mu \rho \upsilon \cdot n\u301b+{\partial}_{i}\left(\u301a2\nu \rho \upsilon \cdot n\u301b{n}_{i}{n}_{j}\right)\text{,}$$ | (53) |

$${\u301a\sigma n\u301b}_{j}={\partial}_{j}\u301a2\nu \rho \upsilon \cdot n\u301b-{\partial}_{i}\left(\u301a2\nu \rho \mathrm{\upsilon}\cdot n\u301b{n}_{i}{n}_{j}\right)\text{.}$$ | (54) |

Further reductions using that κ = ∂_{i}n_{i} and (42) yield the expected jump conditions

$$\u301a\left(\sigma n\right)\cdot n\u301b=-2\u301a\nu \rho \upsilon \cdot n\u301b\kappa ,\text{\hspace{1em}}\u301a\left(\sigma n\right)\cdot T\u301b=2{\partial}_{T}\u301a\nu \rho \upsilon \cdot n\u301b\text{\hspace{1em}}\text{for}\text{\hspace{0.28em}}\text{all}\text{\hspace{1em}}\text{T}\perp n\text{.}$$ | (55) |

Of course they coincide with (18) in the flat case, with (23) in the spherical symmetry, and with (43)–(46) when η = cst (just replace ν with η/ρ). They also coincide with (25)–(26) since **T** is any tangent vector.

## 4 Conclusion

Situations involving a mass transfer across a density jump occur every time a phase change takes place. This occurs in planetary mantles and on the solid-liquid interface of their metallic cores. It occurs also in any situation with phase change or chemical reaction (although in various case a non negligible inertia would have to be taken into account). We have pointed out that a mass transfer across an interface with a density jump is not compatible with a continuous traction vector for a Newtonian fluid. We have illustrated with two simple cases how the new jump conditions appear and shown that the correct conditions lead to a fundamental change in the solutions of the Stokes equations (see Fig. 2). We have derived a jump condition for the traction (25) in terms of a dynamic surface tension, defined in (26) as the opposite of twice the product of the mass transfer $\rho \upsilon \cdot n$ by a coefficient $\u301a\nu \u301b$ depending upon the variation of density and viscosity within the interface. This yields Young–Laplace and Marangoni effects though with a different kind of surface tension.

In a linear elastic solid, we guess that similar arguments would lead to analogous expressions for the jumps, with the displacement instead of the velocity, and the elastic moduli instead of the viscosities. For example seismic propagating and stationary waves could be affected by these jumps at phase change boundaries in the Earth. In the same way, shock waves involve a sharp density discontinuity crossing particules and thus a dynamic surface tension should be considered.

The existence of an infinitesimally thin surface across which a density change occurs is a mathematical idealization. Physically, the transformation takes necessarily place within a finite width. Is it possible to estimate the value of the coefficient $\u301a\nu \u301b$ entering in the equivalent surface tension? It depends upon the evolution of η and ρ within the transition layer. For a univariant phase transition, the interface is a layer of a few molecules depth across which these quantities may not be easy to define, but our approach hypothesizes that the rheology remains viscous and Newtonian. The parameter $\u301a\nu \u301b$ should thus be directly deduced from experiments. For multivariant phase changes like those occurring in a planet (where each silicated phase hosts several type of metallic cations), transitions occur in coexistence regions of 5–30 km width. In this case, the mathematical interface really represents a macroscopic layer of a few ‰ of the planetary radius across which density and viscosity are perfectly defined and could be measured in laboratory experiments. At any rate, the parameter $\u301a\nu \u301b$ cannot a priori be expressed as a function of the density and viscosity values on each side of an interface except in some particular cases. If the density is constant across the interface we have of course $\u301a\nu \u301b=0$. In this case or if the mass transfer is null, the usual continuity of the traction is recovered. If the viscosity is uniform then $\u301a\nu \u301b=\eta \u301a1/\rho \u301b$. We can also show that if the specific volume 1/ρ is an increasing function of q^{3}, and $\eta >0$, then $\mathrm{inf}\left(\eta \right)\u301a1/\rho \u301b\le \u301a\nu \u301b\le \mathrm{sup}\left(\eta \right)\u301a1/\rho \u301b\text{.}$ If η and 1/ρ are linked by a linear relationship then $\u301a\nu \u301b=\frac{1}{2}\left({\eta}^{-}+{\eta}^{+}\right)\u301a1/\rho \u301b\text{.}$ In all these cases, the smaller the jump of specific volume, the smaller the effect of this new jump condition.

Finally, note that if we write the stress field as $\sigma ={\sigma}^{-}{1}^{-}+{\sigma}^{+}{1}^{+}+S\delta \text{,}$ and use eqs. (3), (50), (52), we show that the singular part reads **S** = s**P**_{T} where ${P}_{T}\u2254I-n\otimes n$ denotes the orthogonal projection onto the tangent space to the interface. Thus, **S** is isotropic on the tangent space and verifies the usual properties of a surface tension (Slattery et al., 2007): it is a tangential tensor (i.e. a tensor that maps tangent vectors to tangent vectors) that represents the singular part of the stress tensor and that imposes a traction jump $\u301a\sigma n\u301b=-{\nabla}_{T}\cdot S\text{.}$ The difference between an intrinsic surface tension and the dynamic surface tension **S** is that the latter is proportional to the mass transfer across the interface.

It is beyond the scope of this paper to explore all the consequences of the dynamic surface tension but let us only notice that it might substantially influence inner core convective models results. Indeed the inner core translational convective mode (Alboussière et al., 2010; Monnereau et al., 2010) relies on the fact that **υ** = cst is a solution of Stokes equations together with (**σn**)_{T} = 0 at the boundary. On the contrary, with **υ** = cst the new condition (25) implies

$${\nabla}_{T}\text{\hspace{0.28em}}\left(\u301a\nu \u301b\rho \upsilon \cdot n\right)=0\text{.}$$ | (56) |

At the first order $\u301a\nu \u301b$ and ρ involved in this equation can be considered laterally constant and this condition means that **n** makes a constant angle with the fixed direction **υ**. The surfaces satisfying this condition are called “constant angle surfaces” and this is obviously not the case of a sphere. Thus the translation is not a solution of the Stokes equation in the inner core with a permeable surface involving a dynamic surface tension.

In a forthcoming paper, we intend to generalize the result presented here by relaxing the hypothesis of low Reynolds number. In this case we show that, in the jumps conditions reported above, first we must replace the normal velocity by the relative normal velocity to the interface and second add the classical jump of momentum appearing in the Rankine–Hugoniot conditions.

## Acknowledgments

We thank some colleagues who helped us at an early stage of the work: Élise Poupart, Noé Rabaud and Fabien Dubuffet.