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.

    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;

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

        hpg (pg, q, 0,
    	 ha.x[] += pg);
      delete ({q});

    See also