## Numerical model for bubble plumes

Here we briefly summarize the derivation of our Eulerian-Eulerian model for multiphase (gas and water) flows used to investigate plumes generated by subsurface blowouts in stratified environments. The buoyancy defect at the source is assumed to be due to thermal effects and the presence of a gas phase. The main assumptions are:

1. Boussinesq flow conditions (density variations only accounted in the buoyancy term).
2. Gas bubbles carry a negligible amount of momentum compared to the liquid phase.
3. Constant slip velocity of the gas phase.
4. Linear background stratification.

The starting point for the model derivation, is the mass conservation equation for a phase $k$: $\frac{\partial}{\partial \tilde{t}} \left( \alpha_k \rho_k \right)+ \nabla \cdot \left( \alpha_k \rho_k \mathbf{\tilde{u}}_k \right) =0$

where $\rho_k$, $\alpha_k$ and $\mathbf{\tilde{u}}_k=\left( \tilde{u},\tilde{v},\tilde{w} \right)_k$ are density, volume fraction and Cartesian velocity field for the phase $k$ respectively, $\tilde{t}$ is time and the symbol $\tilde{\phantom{.}}$ denotes dimensional quantities.

Similarly, the momentum conservation can be written as: $\frac{D}{D \tilde{t}}\left( \alpha_k \rho_k \mathbf{\tilde{u}}_k \right)=-\alpha_k\nabla \tilde{p}+\nabla \cdot \alpha_k \widetilde{\boldsymbol{\tau}}_k+\alpha_k \rho_k g \vec{k} - 2\mathbf{\Omega} \times \mathbf{\tilde{u}}_k+\widetilde{\mathbf{M}}_k$

where $\tilde{p}$ is pressure (taken to be equal in each phase ), $\nabla \cdot \alpha_k \widetilde{\boldsymbol{\tau}}_k$ are the molecular and turbulent diffusive terms for phase $k$, $g=-9.8 \, m/s^2$ is the gravity acceleration acting in the vertical direction aligned with coordinate $z$, $\mathbf{\Omega}$ is the rotation vector and $D/Dt=\partial/\partial t +\mathbf{\tilde{u}}_k \cdot \nabla$ is the material derivative.
Momentum transfer between phases is accounted in the $\widetilde{\mathbf{M}}_k$ term and we require: $\sum_k \widetilde{\mathbf{M}}_k=0,$

i.e. phase interactions do not affect the bulk momentum.

As detailed in , under the assumptions presented above, the system of transport equations can be simplified as: $\nabla \cdot \mathbf{u}=0$ $\frac{D \mathbf{u}}{D t} =-\nabla p + \frac{1}{Re} \nabla^2 \mathbf{u} + \left(\theta + Ri \alpha_b\right) \hat{\textbf{k}} - \frac{1}{Ro} \left(\hat{\textbf{k}} \times \mathbf{u}\right)$ $\frac{D \theta}{D t} = \frac{1}{Pe_T} \nabla^2 \theta - \mathbf{u} \cdot \hat{\textbf{k}}$ $\frac{D \alpha_b}{D t} = \frac{1}{Pe_b} \nabla^2 \alpha_b - U_N \frac{\partial \alpha_b}{\partial z}$ $\frac{D \beta}{D t} = \frac{1}{Pe_p} \nabla^2 \beta$

where $\mathbf{u}=(u,v,w)$ is the liquid phase velocity, $\alpha_b$ is the gas volume fraction, $\beta$ is a passive scalar volume fraction, $\theta = T - \zeta z$ is the temperature $T$ perturbation and $\zeta$ is the stratification slope. The system of equations has been non-dimensionalized using the classical velocity, length, time and temperature scales for stratified flows.
Using the inlet buoyancy flux $B_0$ and the buoyancy frequency $N$, these scales are defined as $U_S = (B_0 \, N)^{1/4}$ $L_S = U_S/N$ $t_S = 1/N$ and $T_S = \zeta L_S$ respectively. $U_N = w_s / U_S$ is the non-dimensional slip velocity and $N$ is the buoyancy frequency defined as: $N=\sqrt{-\frac{g}{\rho_0}\frac{\partial \rho_e}{\partial z}} = \sqrt{g \gamma \zeta}$

where $\rho_0$ is the reference density value (taken at $z=0$ for an unperturbed environment), $\rho_e = \rho_0 (1-\gamma \zeta z)$ is the environment density profile and $\gamma$ is the thermal expansion coefficient for the liquid phase.

The non-dimensional groups are the Reynolds number $Re=\frac{U_S L_S}{\nu}$, the Richardson number $Ri=\frac{g \, L_S}{U_S^2}$, the Rossby number $Ro=\frac{U_S}{f \, L_S}$ and the Péclet number for a scalar $m$ $Pe_m =\frac{U_s L_S}{\mathcal{D}_m}$. $\nu$ is the liquid phase kinematic viscosity, $\mathcal{D}_m$ is the diffusivity coefficient for the scalar $m$ and $f=2\Omega$ is the Coriolis frequency (the rotation vector $\mathbf{\Omega}$ is here aligned with $z$).

Detailed derivation of the model can be found in:

 Fabregat Tomàs, A., Dewar, W. K., Özgökmen, T. M., Poje, A. C. and Wienders, N., “Numerical Simulations of Turbulent Thermal, Bubble and Hybrid Plumes”, Ocean Modelling, Volume 90, pages 16-28, 2015. 