# Boussinesq buoyancy

The momentum equation of the multilayer solver becomes \displaystyle \begin{aligned} \partial_t \left( h \mathbf{u} \right)_k + \mathbf{{\nabla}} \cdot \left( h \mathbf{u} \mathbf{u} \right)_k & = - gh_k \mathbf{{\nabla}} (\eta) {\color{blue} - \mathbf{{\nabla}} (h q)_k + \left[ q \mathbf{{\nabla}} z \right]_k} \end{aligned} where the terms in blue have been added and q is the (hydrostatic) pressure deviation due to (small) density variations.

The density variations are described by a scalar field T, which is advected by the flow, and an associated “equation of state” \Delta\rho(T) defined by the user. This gives the additional equations \displaystyle \begin{aligned} \partial_t \left( h T \right)_k + \mathbf{{\nabla}} \cdot \left( h \mathbf{u} T \right)_k & = 0,\\ q(z) & = \int_0^z g \Delta \rho(T) dz \end{aligned}

scalar T;

event defaults (i = 0)
{
T = new scalar[nl];
tracers = list_append (tracers, T);
}

event cleanup (t = end, last)
{
delete ({T});
}

The pressure gradient terms in blue are added to the acceleration of the multilayer solver.

event acceleration (i++)
{

The hydrostatic pressure deviation q is stored on the interfaces between layers (consistently with the Keller box scheme discretisation, see Popinet, 2020). This gives the following vertical discrete integration scheme.

  scalar q = new scalar[nl];

foreach() {
double ph = 0.;
for (point.l = nl - 1; point.l >= 0; point.l--) {
ph += G*drho(T[])*h[];
q[] = ph;
}
}
boundary ({q});

Once the pressure deviation is known, the terms in blue above are added to the face acceleration field ha, using the pressure-gradient macro.

  foreach_face() {
double pg;
hpg (pg, q, 0)
ha.x[] += pg;
end_hpg(0);
}
boundary ((scalar *){ha});

delete ({q});
}