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:

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

The starting point for the model derivation, is the mass conservation equation for a phase :

where , and are density, volume fraction and Cartesian velocity field for the phase respectively, is time and the symbol denotes dimensional quantities.

Similarly, the momentum conservation can be written as:

where is pressure (taken to be equal in each phase [1]), are the molecular and turbulent diffusive terms for phase , is the gravity acceleration acting in the vertical direction aligned with coordinate , is the rotation vector and is the material derivative.

Momentum transfer between phases is accounted in the term and we require:

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

As detailed in [1], under the assumptions presented above, the system of transport equations can be simplified as:

where is the liquid phase velocity, is the gas volume fraction, is a passive scalar volume fraction, is the temperature perturbation and 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 and the buoyancy frequency , these scales are defined as , , and respectively. is the non-dimensional slip velocity and is the buoyancy frequency defined as:

where is the reference density value (taken at for an unperturbed environment), is the environment density profile and is the thermal expansion coefficient for the liquid phase.

The non-dimensional groups are the Reynolds number , the Richardson number , the Rossby number and the Péclet number for a scalar . is the liquid phase kinematic viscosity, is the diffusivity coefficient for the scalar and is the Coriolis frequency (the rotation vector is here aligned with ).

Detailed derivation of the model can be found in:

[1] 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.