Integral models for bubble, droplet, and multiphase plume dynamics in stratification and crossflow

We present the development and validation of a numerical modeling suite for bubble and droplet dynamics of multiphase plumes in the environment. This modeling suite includes real-fluid equations of state, Lagrangian particle tracking, and two different integral plume models: an Eulerian model for a double-plume integral model in quiescent stratification and a Lagrangian integral model for multiphase plumes in stratified crossflows. Here, we report a particle tracking algorithm for dispersed-phase particles within the Lagrangian integral plume model and a comprehensive validation of the Lagrangian plume model for single- and multiphase buoyant jets. The model utilizes literature values for all entrainment and spreading coefficients and has one remaining calibration parameter κ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \kappa $$\end{document}, which reduces the buoyant force of dispersed phase particles as they approach the edge of a Lagrangian plume element, eventually separating from the plume as it bends over in a crossflow. We report the calibrated form κ=[(b-r)/b]4\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \kappa = [(b - r) / b]^4 $$\end{document}, where b is the plume half-width, and r is the distance of a particle from the plume centerline. We apply the validated modeling suite to simulate two test cases of a subsea oil well blowout in a stratification-dominated crossflow. These tests confirm that errors from overlapping plume elements in the Lagrangian integral model during intrusion formation for a weak crossflow are negligible for predicting intrusion depth and the fate of oil droplets in the plume. The Lagrangian integral model has the added advantages of being able to account for entrainment from an arbitrary crossflow, predict the intrusion of small gas bubbles and oil droplets when appropriate, and track the pathways of individual bubbles and droplets after they separate from the main plume or intrusion layer.


Introduction
Jets and plumes occur often in our surroundings, formed due to both manmade and natural reasons. They vary in scale from small, rising hot air plumes formed from smoke stacks to large multiphase plumes created due to volcanic eruptions or subsea oil-well blowouts. Single-phase buoyant jets include sewage water discharge from outfalls (neglecting sediments), hot water plumes from cooling plants, produced water plumes from the petroleum industry, and saline plumes from desalination plants, among others. Multiphase plumes include bubble plumes in lake and harbor aeration, hydrothermal vent plumes laden with precipitating mineral particles, underwater liquid droplet plumes formed in direct ocean CO 2 sequestration, and oil and gas plumes formed due to accidental oil-well blowouts. Jets and plumes attract our attention due to their role in the advection and dispersion of the energy and materials that they carry and the alterations they create in the ambient environment where they are discharged. Due to their small lateral scale compared to the length along their trajectory, buoyant jets and plumes are often modeled using an integral approach based on self-similarity. In multiphase plumes, processes at the particle-scale are further simplified by a discrete particle model for a subset of representative particles. The centerline dilution, plume trajectory, the locations of intrusions, and the pathways of rising particles are common metrics these models need to predict. Since the Deepwater Horizon accident, we have been developing a modeling system for subsea oil and gas plumes based on the discrete particle model [1,2] within an integral plume model framework [3][4][5] that synthesizes modeling approaches in single-and multiphase plumes across the literature [6][7][8]. In this paper, we present this general modeling system for multiphase plumes in stratification and crossflow, its validation to laboratory and field experiments, and we discuss its predictions for canonical test cases of accidental oil-well blowouts in deep water to illustrate the complex role that thermodynamics and mass transfer often plays in multiphase plumes. This work is important to demonstrate the proper applicability range of different modeling approaches from the literature, to highlight some of the limitations of available validation data, and to quantify model performance for predicting important metrics that describe the dynamics and impact of both single-and multiphase plumes in the environment.
In multiphase plumes, an important physics process is that the dispersed phase (bubbles, droplets, or particles-here we follow Clift et al. [9] and use the general term particle) may follow a separate path from the plume centerline due to their own terminal rise (or slip) velocity, and this divergence from self-similarity must be accounted for in an integral model. In Fig. 1, we depict the two main types of dispersed phase separation that occur in multiphase plumes and the formation of secondary plumes and intrusion layers. A plume in pure density stratification (Fig. 1a) rises in the water column, entraining ambient fluid and carrying the dispersed phase upward until a terminal level is reached h P , where the upward momentum of the inner multiphase plume is arrested by the stratification, and the heavy, Environ Fluid Mech entrained water peels, or detrains, from the inner plume, forming a descending outer plume that eventually intrudes, trapping at a level of neutral buoyancy h T [5]. Depending on the slip velocity of the dispersed phase particles, they may either continue through the peeling level, detrain and enter the downdraught plume and later escape, or follow the outer plume into the intrusion layer [10]. Likewise, crossflows, whether in stratified or unstratified environments, may lead to separation of the initial entrained water from the dispersed phase particles (Fig. 1b). Horizontal momentum entrained from the crossflow causes the multiphase plume to bend over in the downstream direction. Separation occurs at the height h S , where the bent plume has been advected farther downstream than the rising particles [11]. After loss of some or all of the dispersed phase, the separated plume of entrained water decelerates in its vertical ascent. In a stratified environment, the separated plume will collapse and form an intrusion (as in Fig. 1b), or in an unstratified ambient with crossflow, the plume transitions to a line thermal. While there are limited data for multiphase plumes with the combined effect of stratification and crossflow, Socolofsky and Adams [11] define plumes as crossflow-dominated when h S \h T and stratification dominated when h S [ h T . Empirical equations predicting these characteristic heights under idealized conditions (i.e., mono-dispersed, non-reacting dispersed phases in linear stratification or pure crossflow) have been fit to the experimental data and were applied by Socolofsky et al. [12] to predict the intrusion height for the Deepwater Horizon accident. However, to include dispersed phase dynamics (compressibility, dissolution, and heat transfer), the combined effects of stratification and crossflow, and to predict the concentrations of dissolved components in the plume and intrusion layers, numerical plume models are required.
Numerical models of the near field of jets and plumes commonly employ an integral approach due their computational efficiency at handling the high length to width ratio (order 10) with adequate accuracy. Only recently have turbulence-resolving computational fluid dynamics (CFD) models become efficient enough to predict multiphase plume dynamics [13,14], and these do not include dispersed phase chemical or thermodynamic transformations. Integral models predict cross-sectionally averaged quantities along the (a) (b) Fig. 1 Schematic diagrams of a a multiphase plume in pure density stratification and b a multiphase plume in stratified crossflow, both showing separation of the dispersed phases from the plume and intrusion formation. In these sketches, h P is the peel height, h T is the trap height, h S is the crossflow separation height, u a is the ambient crossflow velocity, q a ðzÞ is the ambient density stratification and E i ; E a ; E o ; E f , and E s are entrainment fluxes described in Sects. 2

.3 and 2.4)
Environ Fluid Mech plume centerline by applying conservation laws with an entrainment or spreading hypothesis, and in single-phase flows, models have developed in both a Lagrangian [3,6,15,16] and Eulerian [4,17,18] framework. The Lagrangian approach assumes that the plume elements are advected with an average velocity along the plume trajectory, while Eulerian models consider control volumes defined along the plume centerline. These models have evolved to consider three-dimensional trajectories in stratification and crossflow for different release angles and including interactions of multiple jets. Here, we borrow from these models to create the basic numerical framework of our model and to obtain entrainment models that accurately include the different effects of shear entrainment, forced entrainment from the crossflow, and buoyancy effects on entrainment rates. In the multiphase plume literature, two main types of integral models have developed to handle the dispersed phase separation and intrusion formation shown in Fig. 1. The first class of models are the double-plume integral models, which apply to plumes in pure stratification and use an Eulerian framework, solving for an inner, upward-rising plume of dispersed phases and entrained water and a separate, annular outer plume of descending fluid [5,8,19,20]. The inner plume is based on the multiphase model for an unstratified, quiescent ambient by Milgram [21], and the various authors of double-plume models differ in their formulation for exchange fluxes between the inner and outer plumes and the detrainment algorithm that initializes the outer plume; Socolofsky et al. [8] provides a summary and sensitivity analysis of these different approaches. Advantages of the doubleplume model are that it easily lends itself to predict multiple inner plume and intrusion layer sequences when the dispersed phase does not intrude and that it predicts exchange between the inner rising plume and the outer downdraught plume. Disadvantages arise from the fact that it is difficult to allow for intruding particles and that crossflows cannot be included as even very weak crossflows will cause the inner plume to bend in the downstream direction, breaking the symmetry of the outer plume so that it would no longer be annular and thus, making it difficult to formulate the outer plume equations or exchange fluxes between the inner and outer plumes. As a result, the double-plume integral models are strictly appropriate only in quiescent stratification, or where crossflows are weak enough to be ignored, as in the hypolimnia of stratified lakes, where these models were developed to predict weak aeration bubble plumes.
The second class of models solves for the bent plume trajectory in the presence of crossflow using equations similar to the inner plume of the double-plume model. In the single-phase literature, both Eulerian [4] and Lagrangian [3] perspectives are used to predict buoyant jet behavior in crossflows. In the multiphase plume literature, especially for oil-well blowouts, most models are Lagrangian [6,7], and some are based on the spreading hypothesis [22]. In the Lagrangian case, the model solution progresses in monotonic increments along the plume centerline, either bending in the downstream direction or stopping at the height of maximum plume rise h P in the absence of crossflow. Multiphase integral models have evolved to simulate oil plumes in deep water under the action of stratification and crossflow [6,7], to include multiple dispersed phases, including oil and gas [23,24], and to include precipitation reactions at hydrothermal vents [25]. Because of the complex, natural currents in the deep ocean, Zheng and Yapa [23] extended the entrainment hypothesis presented by Lee and Cheung [3] to multi-directional crossflows. Models by Johansen [26] and Chen and Yapa [27] also account for separation of the dispersed phase from the main plume. These models take into account the influence of dispersed-phase particles within the plume and consider them to be separated when they cross the plume edge; however, these authors do not report the numerical approach for tracking the paths of particles within the plume-a complicated task given that the particle and plume dynamics are coupled by the momentum conservation equations. The advantages of these Lagrangian models are their ease at handling realistic crossflows and that dispersed phase particles naturally follow the plume into the intrusion until the dynamics predict that they would separate. Their main disadvantage is that there is no simple means to predict the initial conditions of subsequent plume structures that may form above the first separation point.
Each of these different types of models uses the same set of approaches to account for the dynamics of the dispersed-phase particles. In the recent CFD simulations and the simplest multiphase integral models, dispersed phase particles are inert, having a fixed terminal slip velocity and density, which together predict their net drag (or buoyant force) on a plume element. Fluid particles, especially gases, may change density through compressibility and by dissolution. Most integral models handle these complex, particle-scale dynamics through a discrete particle model [1], which tracks the detailed dynamics for a subset of representative particles and applies their properties to all similar particles in the plume at each computational step. The discrete particle model handles dissolution and stripping [2] by a particle-scale mass transfer equation, and the particle properties, including density and slip velocity, are updated using an equation of state and empirical correlations [9]. Models differ in their complexity, with oil-well blowout models generally requiring non-ideal equations-of-state and solubility models [28]. Here, we utilize a discrete particle model developed for the Deepwater Horizon reservoir fluids [29,30] and focus our discussion on the dynamics of multiphase plumes in the environment.
To compare the different types of multiphase plume models and create a complete modeling suite for oil and gas blowouts, we developed the Texas A&M Oil spill (Outfall) Calculator (TAMOC). The model synthesizes different modeling approaches from singleand multiphase models in the literature and presents a complete discrete particle model together with particle tracking and integral plume simulation modules. In this paper, we apply the TAMOC suite to conduct a thorough validation of the multiphase plume algorithms to data in stratification and crossflow and to compare the results between doubleplume integral models and Lagrangian crossflow models for multiphase plumes in weak crossflows. This paper is organized as follows. In Sect. 2 we present the technical approach to each component of the modeling suite. The methods section also highlights an effective entrainment algorithm for crossflows borrowed from Eulerian single-phase plumes and presents a detailed solution to particle tracking within the Lagrangian integral model in a crossflow. A thorough validation of the double-plume integral model was presented in Socolofsky et al. [8]; here, we validate our crossflow model to laboratory and field data in single-and multiphase plumes in Sect. 3. As an application (Sect. 4), we further compare the double-plume and crossflow integral model approaches using simulations of selected test cases defined by Socolofsky et al. [31] for hypothetical oil-well blowouts in deep water. These simulations demonstrate the role particle-scale thermodynamics and chemistry processes can have on the multiphase plume physics. A summary and statement of conclusions follows in Sect. 5.

Methods
The model presented here is the Texas A&M Oil spill (Outfall) Calculator (TAMOC), which is written in Python and Fortran and freely available under the MIT Public License. The source code for the simulation results in this paper is available from the Gulf of Environ Fluid Mech Mexico Research Initiative Information and Data Cooperative (GRIIDC) [32]. An up-todate version of the model is maintained at Github (http://github.com/socolofs/tamoc). The model is built in an object-oriented way and provides modules for handling the ambient water column data (temperature, pressure, chemical concentrations, and currents), initial particle size distributions, the discrete particle model, and the simulation modules for particle tracking and the integral plume models. TAMOC can be applied to ocean outfalls (referred to as the Texas A&M Outfall Calculator) using the single-phase mode of the model and has been adopted by the U.S. National Oceanic and Atmospheric Administration (NOAA) for near field modeling of subsea oil spills in their General NOAA Operational Modeling Environment in Python (PyGNOME). Simulations are built by utilizing the various methods in each of the main Python modules of TAMOC; the ./ bin/ directory on the Github repository provides example scripts for typical types of model simulation. The model was first introduced in [33,34] and has been applied to simulate the thermodynamic behavior of Deepwater Horizon oil [29,35] and to hindcast a short period of the Deepwater Horizon oil spill [30]. This section presents the technical details of the modeling suite needed to appreciate the validation and application sections that follow.

Discrete particle model: particle-scale dynamics model
The Discrete Particle Model (DPM) computes the physical, thermodynamic, and chemical properties of individual dispersed phase particles. This module is used by all of the multiphase simulation modules in TAMOC. Each dispersed phase particle is defined by its particle type (fluid or solid) and whether or not the particle is reactive or inert. For reactive particles, the particle properties depend on the chemical composition (masses of each chemical component of the particle) and thermodynamic state (temperature, pressure, and salinity in the water surrounding the particle). The accuracy of the methods in the DPM and their implementation within TAMOC have been validated by comparison to measured thermodynamic properties of different gases and petroleum fluids. A rigorous description of the methods in the DPM and a validation to properties of Deepwater Horizon oil is given in Gros et al. [29].
Physical properties of the dispersed phase particles include the equivalent spherical diameter, surface area, particle shape, and slip velocity. These are computed using empirical relations in Clift et al. [9], some of which are summarized in Zheng and Yapa [36]. The thermodynamic and chemical properties include density, viscosity, interfacial tension, heat transfer coefficient, and the fugacity, solubility, diffusivity and mass transfer coefficient of each component of the chemical mixture. The density and fugacity are computed using the Peng-Robinson equation of state with volume translation [37][38][39]. The solubility is predicted by the modified Henry's law equation, in which the partial pressure in Henry's law is replaced by the fugacity, and the Henry's law constant is adjusted from standard conditions to in situ temperature, pressure, and salinity [2,40,41]. Viscosity is estimated from the correlation equation in Pedersen et al. [42], interfacial tension from an equation in Danesh [43], and diffusivities in water from a method in Hayduk and Laudie [44]. Each of these methods require several thermodynamic properties for each chemical component in the mixture, such as molecular weight, critical point pressure, temperature, and molar volume, acentric factor, Henry's law constant, enthalpy of solution, molar volume at infinite dilution in water, Setschenow constant, and the molar volume at the normal boiling point. For simple compounds (e.g., oxygen, methane, benzene), these properties are readily known [45]. For more complex components of petroleum, these Environ Fluid Mech properties can be estimated using group contribution methods, as detailed for our model in Gros et al. [29].
The heat and mass transfer coefficients in TAMOC are expressed as transfer velocities and integrate the fluxes over the whole surface of each particle. Heat and mass transfer rates depend on the shape, rise velocity (wake structure), and internal circulation (convection) of the particle. Fluid particles with internal circulation are called clean particles in the literature because the internal circulations can be shut down by contamination of the particle-water interface by naturally occurring surfactants [9]. These surfactants accumulate on the downstream side of the particle, leading to a concentration gradient of surfactants along the particle-water interface, which results in immobilization of the fluid interface by Marangoni forces. Because clean particles have convection on the inside of the particle, their heat and mass transfer coefficients are generally larger than the transfer coefficients for equivalent dirty particles. The DPM in TAMOC uses mass transfer coefficients for clean and dirty particles from [9,46,47], as described in [35]. Heat transfer coefficients are estimated from the same formulas as the mass transfer coefficients, but with the molecular diffusivity replaced by the thermal conductivity of seawater. These heat and mass transfer functions in TAMOC assume the transfer process is limited by diffusion on the water side of the interface.
The Peng-Robinson equation of state efficiently manages phase transitions. Using this equation of state, the DPM in TAMOC can perform flash calculations to determine gasliquid equilibrium compositions in the two-phase region of thermodynamic state space. Our implementation uses a combination of successive substitution and stability analysis following the method in Michelsen and Mollerup [48]. This allows the model to predict the evolution from live oil (liquid oil with large amounts of dissolved gas) to dead oil (liquid oil with the volatile components removed) via ebullition of the gas out of the liquid phase as a fluid particle decompresses and cools. In general, TAMOC can consider two-phase particles, such as was reported in Gros et al. [30]. Here, we simplify these dynamics and compute the phase equilibrium at the release and initialize only single-phase particles (gas or immiscible liquid) that are assumed to remain single-phase throughout their ascent through the water column (though they may transition from gas to liquid after the volatile components dissolve into seawater). Hence, the DPM in TAMOC is a comprehensive thermodynamic model that provides particle-scale properties at in situ conditions as needed in the simulation models within TAMOC.
The DPM is also designed to be used by other models outside of TAMOC. All of the methods to compute the physical, thermodynamic, and chemical properties are coded as subroutines in a Fortran library. This was done to increase the computational speed, but also to make the algorithms portable to compiled codes. TAMOC provides a Python wrapper to each Fortran subroutine and manages all of the input and output to the Fortran code. In this way, other Lagrangian particle tracking and multiphase models in the literature can adopt the methods in our DPM by interfacing with the Fortran library in TAMOC.

Lagrangian particle model: tracking for reacting bubbles, drops, or particles
The simulation module for particle tracking in TAMOC is the Lagrangian Particle Model (LPM), which is used to track both individual particles released in the water column and the pathways of particles that separate from the integral plume models. The LPM solves for Environ Fluid Mech the advection of a particle along its path (X ¼ xi þ yj þ zk) based on its slip velocity u s k and the ambient currents u a ¼ ui þ vj þ wk according to Here i, j, and k are the unit vectors in the x, y, and z directions, respectively, and u, v, and w are the corresponding ambient velocity components; random displacement resulting from turbulent diffusion is neglected as the model solves for the mean path. The LPM solves for the simultaneous evolution of the particle's chemical composition and thermodynamic state due to mass and heat transfer following the mass and energy conservation laws for the particle given by where subscript i applies to each chemical component in the mixture; m i is the mass, C a;i is the ambient concentration, C s;i is the solubility, b i is the mass transfer coefficient, A is the surface area of the particle, q p is the density of the particle, c p is the specific heat at constant pressure for the mixture inside the particle, b T is the heat transfer coefficient of the particle mixture, T a is the ambient temperature, T p is the particle temperature, DH sol;i is the enthalpy of solution on a molar basis, and M i is the molecular weight. The terms on the right-hand-side of Eq. 3 are the heat transfer between the particle and surrounding fluid and the heat loss from the particle by dissolution, which includes the heat associated with the flux of dissolved mass through the particle interface and the heat of solution of the chemical reaction. In each of these differential equations, t is time, and the initial conditions include the initial particle position X 0 , composition m i;0 , and temperature T p;0 ; q p is computed from m i , T p , and the ambient pressure P using the equation of state in the DPM. When particles contain more than one chemical component, a mass transfer equation (Eq. 2) must be solved separately for each component of the mixture.
The complete system of Eqs. 1 through 3 is coupled due to the dependence of the slip velocity on the particle diameter u s ¼ f ðd e Þ and due to the dependence of the particle diameter on the evolving mixture composition and thermodynamic state. At each time step, the equivalent spherical diameter d e is given from the DPM as and the density q p depends on m i , T p and P. These equations are solved as a coupled system of ordinary differential equations, with all of the instantaneous (u s , A, b i , C s;i , k, q p ) and fixed (c p , DH sol;i , and M i ) particle properties provided by the DPM; ambient data (u a , C a;i , T a , and P) are interpolated via a look-up table in the ambient-properties module of TAMOC. Thus, the LPM simultaneously predicts the fate and transport of multiphase particles as they rise through the water column.

Entrainment models
All of the integral plume models in TAMOC are based on an entrainment assumption, where the flow rate in the plume increases as a result of fluid captured by the plume. Two different types of entrainment are considered. Shear entrainment, or turbulent entrainment, describes fluid engulfed by turbulent structures near the plume edge. This classical type of entrainment is predicted by the entrainment hypothesis, which states that the inflow velocity across the edge of the plume is proportional to a characteristic velocity in the plume [49]. The other type of entrainment is forced entrainment, which results from crossflow intersecting the plume. Forced entrainment is predicted by the intercepted area of the plume and the relative velocity of the crossflow [3].
The proportionality constant in the shear entrainment model is called the entrainment coefficient a s , and different values apply to plumes and buoyant jets and to different definitions of the characteristic plume velocity. For a Gaussian velocity profile with the maximum velocity on the centerline of the plume u c as the characteristic velocity, a general expression for the entrainment coefficient in a single-phase buoyant jet is given by Jirka [4] as where a 1 ¼ 0:055, a 2 ¼ 0:6, / is the angle of the plume centerline from the horizontal plane, and is the local densimetric Froude number, in which b g is the local plume half-width of the Gaussian velocity profile, g 0 c ¼ gðq À q c Þ=q is the reduced gravity of the plume mixture on the centerline, q c is the plume density on the centerline, and g is the acceleration of gravity. The first term on the right-hand-side of Eq. 5 is the pure jet entrainment, and the second term adds the effect of a pure plume.
Equation 5 behaves well in pure jets and plumes in a uniform density ambient, but in stratification, a s is predicted to become infinite as the buoyancy changes sign at a neutral buoyancy point, and F l goes through zero. Jirka [4] proposed a solution to this problem by replacing Eq. 5 for small values of F l by a linear transition region, bounded by F l ¼ AE4:67 (a theoretical plume value), so that the pure plume value a s ¼ 0:083 becomes the maximum allowable entrainment rate. He further extended this to a buoyant jet in crossflow, where / may have any value, and defined the solution for a s in the linear transition region as where a 3 ¼ 0:00131. TAMOC applies Eq. 5 to compute the shear entrainment coefficient when jF 2 l = sin /j [ 21:43 and Eq. 6 otherwise. The value 21.43 comes from the solution of F 2 L = sin / in Eq. 5 when a s takes on the limiting value of 0.083; refer also to Figure 13 in Jirka [4]. Because the integral models in TAMOC use a top-hat velocity profile, a s predicted by Eqs. 5 and 6 is multiplied by ffiffi ffi 2 p to convert the Gaussian entrainment coefficients to equivalent coefficients for a top-hat model. Similar conversions are required to compute F l from the top-hat variables.
The forced entrainment in TAMOC is computed from the crossflow velocity and plume geometry and represents the total amount of ambient fluid intercepted on the windward side of the plume. Lee and Cheung [3] present the equations for forced entrainment for an Environ Fluid Mech arbitrary control volume. The crossflow also breaks down the symmetry of the shear entrainment as the plume bends over and forms a vortex pair. In this region, the total entrainment rate in a crossflow is usually taken as the greater of either the shear or forced entrainment rates [3]. For different definitions of the shear entrainment, the transition from shear-dominant to crossflow-dominant entrainment can cause model instabilities. Different solutions to this problem are presented in Lee and Chu [16] and Lee et al. [50]. We found that the solution for shear entrainment by Jirka [4] in Eq. 6 to handle F l ! 0 also solves this problem, yielding a smooth transition from shear-to crossflow-dominated entrainment. Hence, TAMOC always computes the shear entrainment from Jirka [4] and the forced entrainment following Lee and Cheung [3], taking the maximum from these two entrainment types as the local total entrainment.
The equations in Lee and Cheung [3] for the forced entrainment assume that the crossflow is in the x-direction. In general, the crossflow can be in any direction, and in TAMOC, we allow arbitrary currents, with velocity components u a , v a , and w a in the x-, y-, and z-directions, respectively, where z is in the fixed reference frame of the discharge and is directed positive downward. To use the solution in Lee and Cheung [3], we apply a local coordinate rotation. Let h 0 be the horizontal angle from the crossflow direction and / 0 the vertical angle from the crossflow direction. These are related to the TAMOC coordinate system by h 0 ¼ h À h a and / 0 ¼ / À / a , where h is the angle measured from the x-axis in TAMOC, / is the angle measured from the horizontal plane in TAMOC, and h a and / a are given by Hence, the equations in Lee and Cheung [3] for a unidirectional crossflow can be used for any arbitrary crossflow by the above transformation, applied locally at any plume control volume. Finally, all of the entrainment models described in this section were derived for singlephase buoyant jets. Socolofsky et al. [8] calibrated a double-plume integral model for a multiphase plume and found that the best-fit entrainment coefficient for the inner plume (a i ) was 0.055, about half of the single-phase value for the top-hat model 0:083 ffiffi ffi 2 p ¼ 0:12. Milgram [21] reported measured entrainment coefficient values for a wide range of bubble plume experiments, with values spanning 0.057-0.23 for top-hat models, the results depending on a bubble Froude number. Milgram [21] specifically compared the measured data to correlations in the form of Eq. 5 and concluded that the values in bubbles plume did not follow the same model. Hence, the true value of the shear entrainment rate for multiphase plumes is quite complicated, not expected to be equal to that of single-phase plumes, and is unknown for multiphase plumes in crossflows. On the other hand, the values of the entrainment coefficient are the same order of magnitude as those predicted by singlephase equations (about 0.1). To reduce the number of free parameters in TAMOC, we assume the entrainment models for multiphase plumes can be given by the single-phase equations presented here. The fact that this approximation yields acceptable results is borne out in the Validation, Sect. 3.

Stratified plume model: multiphase plume model for quiescent or near still ambient environments
The module for the Stratified Plume Model (SPM) follows the equations presented in Socolofsky et al. [8] and uses an Eulerian framework to predict the steady-state solution for multiphase plumes in stratified and still, or nearstill, environments. This model is similar to Asaeda and Imberger [5] and McDougall [19], but with a continuous peeling equation, as introduced by Crounse et al. [20]. The SPM in TAMOC extends the model in Socolofsky et al. [8] by providing a non-ideal equation of state in the DPM, by applying Eqs. 5 and 6 to compute the entrainment coefficient for the inner plume, and by allowing more than one simultaneous particle size or type within the plume. When multiple particles are simulated, the total buoyant force is computed by summing over the contribution from all particles present. The detrainment algorithm in the SPM is continuous in the sense that the peeling flux from the inner plume into the outer plume takes on different values at each height along the inner plume. Yang et al. [13] proposed a new formulation to compute the peeling flux E p based on post-processing of the velocity field computed from a large eddy simulation (LES) model of bubble plumes in quiescent stratification, which also results in continuous peeling behavior-hence, this is a real, physical phenomenon. If particles are allowed to move from the inner plume to the outer plume with E p , then the computational efficiency of the DPM will breakdown for the outer plume: at each inner plume grid point, a new set of dispersed phase particles would be initialized in the outer plume. This is further complicated by the fact that no algorithm exists to determine when an inner plume particle would be transported into the outer plume. Due to the present lack of process understanding and available validation data, the SPM assumes that dispersed phase particles remain in the inner plume and do not peel into the outer plume.
Since particles in the SPM remain in the inner plume, if an intrusion forms before the plume reaches the surface, multiple inner and outer plume intrusion structures will form. For most gas bubbles, it is a good assumption that the gas remains in the inner plume, yielding good validation to measured laboratory and field data [8,51]. Particles having small slip velocities (e.g., very small particles, or particles that are nearly neutrally buoyant) may tend to intrude with the peeling fluid [10]. When using the SPM, these are best introduced as part of the continuous phase. To track these particles into the intrusion layer, a different type of model is needed. In the case of an ambient crossflow, the bent plume model can provide this capability.

Bent plume model: multiphase plume model for stratified crossflow
To simulate a multiphase plume in stratified crossflow, the Bent Plume Model (BPM) module of TAMOC adapts the methods in Lee and Cheung [3], Lee and Chu [16], and Jirka [4] using a Lagrangian framework to solve for the steady-state solution and adds the dispersed phase following the approaches in Johansen [7,26] and Chen and Yapa [52]. Some of these approaches were also used by [53,54] to simulate blowouts like the Deepwater Horizon. When dispersed-phase particles are added, their buoyancy must be added to the net plume buoyancy, and there needs to be an algorithm to track the particles within the plume, allowing them to separate from the entrained plume fluid under the actions of stratification and crossflow (see Fig. 1b).
Environ Fluid Mech

Momentum conservation for multiphase plume
When particles are introduced to the BPM, their properties are provided from the DPM, and the net buoyant force on the plume resulting from the dispersed phases must be included in the conservation of momentum. For a single-phase Lagrangian element, the change in vertical momentum d(Mw) / dt is equal to the net vertical buoyant force on the Lagrangian element F b , given by where M and w are the mass and vertical velocity of the Lagrangian plume element, q r is a reference density, q a is the density of the ambient fluid at the present location of the Lagrangian element, and q is the density of the entrained fluid in the Lagrangian element; in TAMOC, we take the vertical coordinate to be depth, which is positive down, and we set q r to the density in the middle of the water column. When a dispersed phase is introduced, we assume that the void fraction is very small, so that the momentum of the dispersed phase itself is negligible. This is the normal assumption for a bubble plume and is valid as long as the source conditions are plume like and do not form a multiphase jet. Then, the buoyant force of the dispersed phase F p acting on the Lagrangian element results from the mass of the displaced volume and reduced gravity of the particle and is given by where c is the momentum amplification factor accounting for turbulent kinetic energy production by the dispersed phase [21], j is a reduction factor between 0 and 1 based on the particle position relative to the plume centerline (see next subsection), M p;i is the total mass of particles of type i in the Lagrangian element, and M p;i =q p;i is the volume occupied by particles of type i. Thus, qM p;i =q p;i is the mass of the fluid in the plume displaced by the dispersed phase particles. When we include entrainment of ambient fluid _ m e having a vertical velocity w a , the complete conservation of vertical momentum equation becomes The remaining elements of the BPM include conservation of mass, horizontal momentum, salt, heat, dissolved component masses, and passive tracers. Mass and heat exchange between each dispersed phase particle and the plume fluid are predicted by Eqs. 1-3; the other conservation equations are solved in a standard Lagrangian plume model approach [3,16].

Tracking dispersed-phase particles in a multiphase plume
Because dispersed phase particles have a slip velocity relative to their surrounding fluid, when a plume bends over in a crossflow, particles in the plume will not follow the same trajectory as the plume centerline. Eventually, particles may separate from the plume on the upwind side (see Fig. 1b), which will remove buoyancy from the plume. Lagrangian plume models developed for oil-well blowouts have been adapted to track gas bubbles and oil droplets by Chen and Yapa [52] and Johansen [26], and their algorithms for particle Environ Fluid Mech separation have been validated to experimental data in Socolofsky and Adams [11]. Johansen [26] does not describe the details of the particle tracking, and Chen and Yapa [52] only explain that the entrainment velocity toward the plume centerline should be included for the advection of the particles in the tracking equations. Because the vertical momentum conservation (Eq. 11) depends on whether the particles are in the plume or not, particle tracking must be conducted simultaneously with the plume solution, and the algorithm for particle tracking is not trivial. For these reasons, we provide the details of the particle tracking algorithm in TAMOC. Particle tracking in the BPM involves the simultaneous solution of multiple Lagrangian elements (particles and the plume element), each with different advection velocities. Due to their differential advection, particles and the plume element will travel different distances in a fixed numerical time step. In Fig. 2 we depict the Lagrangian plume element at time t 0 and after one solution time step t 0 þ Dt along with the trajectory of one dispersed phase particle. The figure also defines the local coordinate system at the base of a plume element ðk; n; gÞ, with unit vector l along the plume centerline line, m pointing at the center of curvature, and n normal to l and m. The global coordinate system has unit vector s along the plume centerline and / the angle from the horizontal. A fixed Cartesian coordinate system (x, y, z) is also defined with origin at the water surface directly above the discharge point (z is positive down); U is the local speed of the plume element in the s-direction, À u s k is the vertical slip velocity of the dispersed phase particle in the Cartesian coordinate system, and u e is the velocity resulting from the entrainment, directed toward the plume centerline. The local half-width and height of the plume element are b and h.
If we denote _ m e as the local mass inflow rate of entrainment and assume the entrainment velocity decreases linearly from the value at the plume edge to zero on the plume centerline, then in the local coordinate system Our assumption of a linear profile for entrainment rate is approximate. The true lateral velocity profile for a plume in quiescent conditions can be obtained from continuity and is zero at the plume centerline, directed outward and increasing for small radius, and directed inward and decreasing toward the plume edge. In crossflow, the radial symmetry breaks down, with the formation of a vortex pair on the downstream edge of the plume, and when forced entrainment is dominant, the entrainment profile is further complicated by the capture of upstream fluid. The main effect of the entrainment velocity used here is to delay the expulsion of particles from the plume due to the inward sense of the entrainment velocity vector, and we assume this delay is caused by turbulent eddies in the plume that trap the particles for some time. Because of the complicated nature of the true lateral velocity profile and our desire to have an analytical advection solution, we propose the linear profile used here. The validity of this assumption is tested later in the validation sections.
The vertical slip velocity of the dispersed phase particle can also be projected onto the local coordinate system, defining If we further define the prefactor to the entrainment velocity as f e ¼ _ m e =ð2pb 2 hq a Þ, then the equation for the dispersed phase advection in the local coordinate system of a plume element is where X ¼ kl þ nm þ gn. Because the dispersed phase particles will move a different distance in a time Dt than the Lagrangian element (e.g., they may pass completely through the Lagrangian element and into neighboring elements; see Fig. 2), this advection equation cannot be coupled to the BPM conservation equations using a single numerical time step.
We solve the problem of these two advection time scales by solving Eq. 14 analytically at each time step and then determining where the particle trajectory intersects the base of the Lagrangian element at the new position at time t 0 þ Dt. During each model time step, the local variables U, u s , and f e are taken as constant. The initial conditions are ðn 0 ; g 0 ; k 0 Þ, the position of the particle at the base of the Lagrangian element at time t 0 . Then, the analytical solution to Eq. 14 is In the numerical solution, we solve for the time t 0 ¼ t 0 þ dt when the particle intersects the base of the new Lagrangian element at time t þ Dt. This usually shorter time step dt is then used as the particle time scale in all particle equations (including mass and heat transfer). This approach works because the Lagrangian plume models are solving for the steady-state solution.
Particles are tracked with the plume until their radial position r from the plume centerline is greater than the half-width b. After particles exit the plume, they can be tracked using the LPM, but no new plumes are assumed to form (as, e.g., in the SPM, where multiple intrusion layers may form). As particles move toward the edge of the plume, it is also reasonable to reduce their effectiveness at contributing buoyancy to the plume, as was done by Lai et al. [55] for sediment clouds. This is because particle separation in real plumes is gradual (note that Socolofsky and Adams [11] find a separation point by interpolating between points before and after separation), and this helps to avoid a discontinuity in F p when particles leave the plume. The coefficient j in Eq. 10 provides the buoyancy efficiency, and we considered several functional forms in its calibration (see Sect. 3.2 below). The best functional shape we found for j was Hence, particle tracking in TAMOC utilizes analytical solutions to the advection equation and gradually reduces the influence of particle buoyancy inside the plume until particles exit through the plume edge.

Initial conditions
The three models in the TAMOC suite that solve coupled sets of ordinary differential equations (ODE) are the LPM for tracking individual bubbles, drops, or particles and the two plume models, the SPM and BPM. For particle tracking, the initial conditions are the particle position, its composition on a mass-per-component basis, and its temperature. All properties are computed from the DPM, and ambient properties of the water column are handled by an ambient data module. For the plume models, the initial conditions include the fluxes, composition, sizes, and temperature of the dispersed phase particles and the release position, orientation, and diameter (all releases are treated as circular). For releases that do not include any ambient fluid, there needs to be an estimate of the initial amount of ambient water that mixes with the plume in the zone-of-flow-establishment; otherwise, the initial volume and momentum fluxes of the entrained fluid are undefined, and the solution fails. In these cases, we use the bubble plume Froude number initial condition proposed by Wüest et al. [56] and detailed for the SPM in Socolofsky et al. [8]. Alternatively, the zoneof-flow establishment could be considered, as in Premathilake et al. [57]. If the release is single-phase or if the release includes ambient fluid as part of the release, then the velocity of the released fluids substitutes for the Wüest et al. [56] initial condition.

Numerical solution and model framework
The model equations for each of the simulation modules (LPM, SPM, and BPM) are stiff, and we solve each set of equations using the same numerical techniques. The LPM equations are stiff due to discontinuities in some of the bubble properties (e.g., slip velocity has a discontinuity when the bubble shape changes from spherical cap to elliptical and to spherical) and for multi-component bubbles due to the rapid composition changes that occur as one of the components approaches very low mass fraction. The two plume models, the SPM and BPM, are stiff because changes are very rapid at the release compared to the main body of the plume and because state variables undergo abrupt changes at neutral density points caused by the stratification. Each of these models is coded in Python, and we use the SciPy ODE integrator based on the Variable-coefficient ODE (VODE) solver with backward differentiation formulas (BDF) for stiff problems (documented at https://docs.scipy.org/doc/scipy/reference/ generated/scipy.integrate.ode.html as of March 22, 2018). This SciPy module is a Python port of the Fortran VODE solver in the Netlib repository (DVODE). We use a 5thorder integrator with an adaptive step size. Because of some of the discontinuities in the property data for particles, we had to adjust the normal error criteria for the adaptive step-Environ Fluid Mech size solver; we use an absolute tolerance parameter (atol) of 1 Á 10 À6 and a relative tolerance parameter (rtol) of 1 Á 10 À3 . These values produce mostly stable solutions with a balance between the number of calculation steps and an acceptable level of numerical accuracy.

Model validation
The SPM in TAMOC solves the same equations as in Socolofsky et al. [8], and they present a comprehensive calibration and validation of that model. We have tested the model implementation in TAMOC to ensure that the SPM reproduces the results in Socolofsky et al. [8], and this is not repeated here. The equations of state in the DPM are validated comparing to data for the Deepwater Horizon oil in Gros et al. [29], and the LPM has been validated to analytical solutions for particle advection. The BPM remains to be validated, and we present the results of the validation in this section with a focus on the entrainment models, particle tracking, calibration of the buoyancy parameter j, and the numerical solution using an adaptive time step, stiff equation solver.

Single-phase plume validation
To validate the single-phase behavior of the BPM, we primarily use the validation cases for CorJet reported in Jirka [4]. All of the measured data used in these validations are available for download from the CORMIX website at http://www.mixzon.com/benchmark/ as of March 22, 2017. The JetLag model by Lee and Cheung [3] has also been compared to this validation database in Lee and Cheung [3], Lee and Chu [16] and Lee et al. [50]. Overall, CorJet and JetLag yield similar performance, with each model having its own advantages. The BPM shows a similar level of performance since it is based on the same modeling approach as JetLag, but with the shear entrainment algorithms of CorJet. The main differences in performance for the BPM relative to these other two models stem from the shear entrainment formulation and the adaptive step size numerical solver; both CorJet and JetLag have used constant step size integration methods in the cited references. We show a few cases of these validation exercises here, each selected to highlight the subtle differences present in the BPM solution.
Since CorJet uses Gaussian profiles, the validation database above presents measured data in a form appropriate for comparison to a Gaussian profile. The Lagrangian-based models, such as the BPM, solve for a tophat profile. We use standard conversions between these models throughout this section (e.g., Fischer et al. [58]) to present the results from the BPM as those for equivalent Gaussian profiles. In this way, the figures presented below can be directly compared with figures in Jirka [4].
The simplest case is a vertical buoyant jet into a still, uniform-density environment. When a single-phase source enters the receiving water, it carries a finite amount of momentum, initially following asymptotic relationships for a jet. After a distance proportional to the jet-to-plume transition length scale L M ¼ M 3=4 0 =B 1=2 0 , the discharge behaves like a plume, where M 0 is the initial kinematic momentum flux (u 2 0 pD 2 =4) and B 0 is the initial kinematic buoyancy flux (u 0 ðq a À q 0 ÞgpD 2 =ð4q a Þ). D is the diameter of the discharge port, the subscript 0 denotes cross-sectionally averaged values at the discharge location, and the subscript c in the figures and equations below denotes values of a Gaussian profile on the buoyant jet centerline. In Fig. 3 we depict the non-dimensional centerline velocity as a function of non-dimensional distance z=L M from the discharge for a vertical buoyant jet in a stratified environment. The corresponding centerline dilution S c is plotted in non-dimensional space in Fig. 4; here, F 0 ¼ u 0 = ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðq a À qÞ 0 gD=q a p is the densimetric Froude number at the discharge. From Fig. 3 we conclude that the BPM transitions from a jet [momentum conserving with constant M . This behavior means that the model slightly over-predicts the centerline velocity in the plume region, and this performance can be compared with Figure 7 in Jirka [4], in which case the CorJet model shows earlier and more significant under-prediction of the non-dimensional momentum flux. The centerline dilution in Fig. 4 also indicates that the asymptotic behavior of the model in the jet and plume regions as well as the jet-to-plume transition is captured correctly. The BPM performance for centerline dilution is nearly identical to CorJet, owing to the fact that we use the same shear entrainment algorithm. The over-prediction of the centerline velocity of the plume can be corrected by changing the entrainment coefficients; however, plume models generally give greater importance to dilution than momentum, which explains why the present entrainment model, which matches better the centerline dilution than centerline velocity, is considered optimal.
To further test the entrainment model and to validate the plume tracking equations, we evaluate the trajectory of buoyant jets discharged at an angle into still environments. In this case, the experiments used dye visualization, and the outer edge of the dyed region was digitized for comparison to the model results for centerline trajectory and plume halfwidth. We compare the model trajectory to that of a desalinization brine discharge in a uniform ambient in Fig. 5. The plume is discharged at 55 to the horizontal. The negatively-buoyant discharge decelerates away from the outlet and eventually changes direction, accelerating and falling back to the level of the discharge. This test validates the linear transition of the shear entrainment between decelerating and accelerating buoyant jets given by Eq. 6. The BPM performs similarly to CorJet, with better agreement for x=D [ 80. We attribute this to the adaptive-step and stiff equation solver, which better preserves the initial momentum near the source, where a very small step size is used.
We test the BPM performance in crossflows initially in an unstratified environment. In Fig. 6 we report the centerline trajectory of a pure jet discharged at different angles into a uniform crossflow. In this figure, / is the angle from the horizontal and h is the angle relative to the crossflow in the horizontal plane. U a is the ambient crossflow velocity and U j is u 0 . The centerline trajectory of jets discharging in the crossflow direction (h ¼ 0 ) have excellent correspondence with the measured data, and those discharged against the crossflow (h ¼ 180 ) show degraded performance after turning into the crossflow direction, but do capture well the maximum upstream penetration of the jet. The model performance is at a similar skill to CorJet, but differs due to the different forced entrainment used in the BPM, which follows Lee and Cheung [3] instead of Jirka [4]. This different model for the crossflow entrainment allows the BPM to more closely match the data farthest from the source compared to CorJet (Figure 23 in Jirka [4]). The degraded performance for all models for jets in opposition to the crossflow results from the loss of selfsimilarity for real jets in these cases; hence, the current performance level is remarkable. Finally, we compare the BPM performance for centerline trajectory in stratified crossflow. For this evaluation, we use data from Hirst [17], which was not used by Jirka [4], but has been used for the Clarkson Deep Oil and Gas (CDOG) model in single-phase mode, reported in Zheng and Yapa [64]. In Fig. 7, we show the model centerline trajectory for two different buoyant jets released vertically into a stratified crossflow. In both cases, the BPM under-predicts the vertical rise rate of the buoyant jet. The BPM rises slower than results reported for CDOG in Zheng and Yapa [64] and also lies below predictions for simulations we conducted using the JetLag model; both CDOG and JetLag are a closer match to the data than TAMOC. We attribute these differences to our adaptive step-size and stiff equation solver, which takes very small steps close to the discharge, likely entraining more ambient fluid from the dense water near the source than in the CDOG and JetLag models. Whereas, this better numerical solution close to the source benefitted TAMOC in a uniform ambient (i.e., results in Fig. 5), in stratification, our solution may lead to over-prediction of the early stage entrainment. Nonetheless, this performance Environ Fluid Mech remains close to that of other models, and these validations confirm that the BPM correctly solves the buoyant jet equations and achieves an averaged performance on par with similar single-phase models in the literature.

Multiphase plume validation: laboratory experiments
Many papers in the literature report laboratory observations for pure bubble plumes (releases of air bubbles from a straight orifice, a nozzle, or porous plate), and validation of the SPM to data in a quiescent ambient is presented in Socolofsky et al. [8]. Fewer studies consider simultaneous discharge of multiphase fluids, and here we validate the BPM to data for bubble plumes and bubbly and oily jets in crossflow environments.
Socolofsky and Adams [11] present an important dataset for pure bubble plumes in crossflow which was also used by Chen and Yapa [52] and Johansen [26] to validate their bubble separation algorithms. Socolofsky and Adams [11] present data for the trajectory of the entrained ambient water and gas bubbles and use these data to correlate the separation height h S with the crossflow and bubble plume parameters. In Fig. 8, we present a typical comparison of the BPM simulation result with the trajectory data for a bubble plume in a modest crossflow. The BPM solution for the entrained plume fluid matches well with the digitized edge of the dye visualization from the experiments. The trajectory of the bubbles is slightly ahead of the region over which the bubble column was observed in the laboratory; yet, the separation height is very similar between the model and observed data.
We compare the simulated separation heights for all of the bubble plume experiments in Socolofsky and Adams [11] in non-dimensional space in Fig. 9. Socolofsky and Adams [11] report that the observed separation height is somewhat ambiguous due to the turbulent boundary of the injected dye tracer and finite width of the bubble column. Hence, in the figure we connect the range between the measured minimum and maximum separation heights reported in their paper by lines terminated by a ? symbol. To give the best overall match between the model and the data, we adjusted the power-law relationship of j, the parameter that reduces the buoyancy of the dispersed phase as it approaches the plume edge. The form for j reported in Eq. 18 was the best-fit that produced the model results in the figure. To extend the scope of the validation runs, we compared the model simulations to the measured data for runs with the same parameters as in the experiments and for selected simulation runs with field-scale initial conditions in the model, but plotted in the same nondimensional space. For separation heights with u 1 =ðB=zÞ 1=3 [ 0:5 or u s =ðB=zÞ 1=3 \3 (modest bubble plumes with modest to low separation heights), the model predicts values that fall within the measured range for most experiments. For weaker bubble plumes that separate far from the source, the model is more likely to predict u s =ðB=zÞ 1=3 to be greater than the measured range: the model over-predicts the separation height in these cases. However, the separation height is given by the moment that the bubbles cross the halfwidth boundary b of the modeled plume. Because j decays as ððb À rÞ=bÞ 4 , the buoyant force of the bubbles is substantially removed from the plume well before the final separation height, and the model continues to predict the trajectory of the entrained fluid well, similarly to that shown in Fig. 8.
A similar type of experiment for bubbly jets is reported in Zhang and Zhu [65]. Their experiments include air and water premixed and discharged vertically from a 6 mm diameter nozzle into a uniform crossflow. We compare BPM simulation results to observations for four of their experiments in Fig. 10. For a multiphase discharge, there is some ambiguity regarding the initial velocity (hence, initial momentum flux) of the discharged water owing to the unknown void fraction at the nozzle exit. Zhang and Zhu [65] considered three different formulations for the initial water jet velocity, and the results in Fig. 10 use the superficial water velocity u 0 ¼ Q=ðpD 2 =4Þ, where Q is the flow rate of water, as the initial condition. For these experiments, the model results are weakly dependent on the water exit velocity within the range of values reported for each experiment by Zhang and Zhu [65], and our optimum results are obtained for this lowest possible initial jet velocity, likely due to entrainment of ambient water that occurs in the zone-offlow-establishment, giving an effective initial jet velocity apparently close to the water superficial velocity at the nozzle. In each case, the trajectories of the bubbles outside the water jet overlap with the measured data, showing slightly steeper slopes. These simulations only include the mean bubble size as reported in [65], and other slopes would be obtained for different bubble sizes for the whole bubble size distribution. Fig. 9 Non-dimensional plot of crossflow separation height h S with crossflow velocity u 1 compared to data: dots are modeled separation heights, dashed and dotted lines are best fits reported in Socolofsky and Adams [11] to the laboratory data, and the plusses connected by lines show the measured range from the minimum to the maximum separation heights reported in Table 1 in Socolofsky and Adams [11] Environ Fluid Mech Likewise, the centerline of the water jet is in good agreement with the observed data, but with the BPM consistently under-predicting the rise rate of the jet far from the bubble column. The under-prediction of the BPM is greatest for the weakest crossflows relative to the bubbly jet buoyancy flux. Our observations for pure bubble plumes confirm that the model is likely to under-predict the separated jet rise rate in these cases. This results from the fact that water that is advected through the bubble column above the separation height (e.g., at z ¼ 0:15 m in experiment B33 in Fig. 10) receives some vertical momentum from the bubble drag; hence, there is a non-zero vertical velocity in the wake of the bubble column above the separated plume (i.e., in the region approximately bounded by x [ 0:3 m and z\0:2 m in experiment B33 in Fig. 10). The integral model does not have information about this background, vertical velocity field; thus, it will under-predict the vertical rise of the far-field separated jet. In strong crossflows, there is less time for the bubble column to impart vertical momentum to the wake, and this effect is much less, as observed for the much better model agreement in experiments C33 and C35 in Fig. 10.
Quantitative observations for the trajectory of bubble plumes or bubbly jets in stratified crossflow have not been reported in the literature. To fill this gap, we conducted a set of dye visualization experiments similar to the methods in Socolofsky and Adams [11]. In Fig. 11, we present a comparison of the BPM for one experiment in stratified crossflow. The stratification profile is described by the buoyancy frequency N equal to where q is a reference density and z is the depth (positive down). Using the correlation equations in Socolofsky and Adams [11,51], the separation point for this plume is predicted to be 0.28 m depth, and the peel point in quiescent stratification would be above the free surface. Although Socolofsky and Adams [11] suggested h S \h T as the criterion for a plume to be crossflow dominated, our recent experiments and the paper by Socolofsky et al. [12] prefer the criteria h S \h P , and we follow this definition here. Hence, the plume in Fig. 11 is crossflow-dominated. The results in the figure demonstrate that the simulated trajectories of the bubble column and the entrained plume fluid agree with the trends in the observations. The separation height for the simulation is similar to the empirical value (0.23 m depth in the simulation). After the entrained fluid separates from the particles, it descends due to its excess buoyancy to form an intrusion in the stratified background current near 0.4 m depth. The simulated intrusion dynamics and trap depth are in good correspondence with the area outlined by the dye in the experiments, though it is difficult to judge the bottom edge of the intrusion boundary in the experiments due to the high dye dilution in this region.
As a final laboratory-scale validation case, we consider data for an oily jet reported in Murphy et al. [66]. In their experiments, the release orifice was cosine-shaped, leading to a contraction of the flow shortly after emission from the orifice. We calibrated a contraction coefficient for this nozzle by comparison to their results for a single-phase jet. The results shown in Fig. 12 were obtained with a contraction coefficient of 0.4; the single-phase jet is shown in panel d. These experiments were designed to test the effects of oil breakup with and without dispersant; hence, the oil velocity at the orifice was large. Because the released fluid is pure oil, all of the momentum creating the jet comes from the dispersed oil phase. The default behavior of the model, as described above in Sect. 2.5, is to neglect the momentum flux of the dispersed phase. This assumption is valid as long as the jet-to-plume transition length scale l M of the dispersed phase is small. In the oily jet experiments of Murphy et al. [66], l M given by the oil momentum and buoyancy flux is not negligible, and we modified the BPM to include the momentum of the discharged oil. This was not Fig. 11 Trajectory of an air bubble plume with an airflow rate of 0.6 Nl/min and bubble slip velocity of 0.17 m/s into a stratified crossflow (buoyancy frequency N ¼ 0:61 s À1 and crossflow velocity u a ¼ 4 cm/s) Environ Fluid Mech required in the Zhang and Zhu [65] experiments because the momentum in their bubbly jets resulted from the water in the release, which is always included in the TAMOC plume models. The resulting simulations shown in Fig. 12 include the momentum of the oil phase and show good agreement for the centerline trajectory of the oily jets and the single-phase jet.
The trajectories of oil droplets in the jets are also well matched. In panel a, the oil discharge is without chemical dispersant injection, and the oil droplets are largest. In the figure we annotate the trajectories for oil droplets with initial diameters of 5, 2.3, and 1.1 mm. The inset showing the sizes of oil droplets near the 2.3 mm trajectory line is from the experiments and shows that similar sized droplets are found in the observations there (in the 2-4 mm range). In panels b and c, chemical dispersants are injected, giving rise to smaller oil droplets that do not leave the buoyant jet; TAMOC also predicts that the oil remains inside the jet in agreement with the observations. These experiments confirm that the algorithms for particle tracking in the BPM yield reasonable results for predicting the conditions under which oil droplets will rise out of a plume or stay inside the plume.
These validation cases for laboratory experiments of bubble plumes and bubbly and oily jets demonstrate that the bubble tracking algorithms in the BPM are working correctly and that the buoyancy reduction factor j gives good agreement for both the separation height and the trajectories of the bubbles and entrained plume fluid in bubble plumes and bubbly and oily jets in crossflow. The limited data for stratified crossflow also show that the BPM is an appropriate tool in crossflow-dominated conditions, yielding good agreement to the measured data for the maximum rise height of entrained plume fluid and the intrusion depth.

Multiphase plume validation: feld experiments
An important validation dataset for field-scale multiphase plumes is the DeepSpill experiment, conducted by Sintef, Norway, in the Norwegian Sea. The experiment simulated an accidental oil-well blowout and involved the release of natural gas (approximately 99% methane) together with seawater, marine diesel, or crude oil from a multiphase discharge at 844 m depth [67]. The primary subsurface observations were of the oil droplet and bubble size distribution at the release using cameras on a remotely operated vehicle (ROV) and of the distribution of gas bubbles in the water column from acoustic backscatter using the echo sounder on the support vessel. Integral models by Johansen [26] and by Chen and Yapa [27] have been previously validated to this classical dataset.
Here, we compare the BPM simulation results for the trajectories of bubbles to the measured acoustic data for two experiments involving hydrocarbon releases. In Fig. 13, we present results for a release rate of 0:7 Nm 3 =s of natural gas and 60 m 3 =h of seawater, which formed a bubbly jet of gas and water. The colored blocks in the figure are the acoustic backscatter intensities measured from the echo sounder [67,68]. The crossflow data are taken from a moored acoustic Doppler current profiler (ADCP) that was recording during the experiments. The simulation results show that the entrained fluid separates from the bubbles near the source (maximum rise depth of about 700 m), and descends to form an intrusion in the background stratified crossflow. Several bubble sizes and their dissolution were simulated from the measured bubble size distribution, and their trajectories are given by the dash-dot lines. These lines correspond nicely with the central region containing bubbles identified by the echo sounder, and the maximum rise height of the bubbles, where they are considered fully dissolved in the model (about 200 m depth), also agrees with the observations.
In Fig. 14, we present simulations for the experiment involving natural gas discharged with crude oil. The currents were similar during the crude oil experiment compared to the natural gas and seawater experiment, and the plume is crossflow dominated, trapping close Fig. 13 Model comparison with the DeepSpill experimental data for natural gas release (dashed lines are bubble paths; continuous lines are plume edges and centerline). Orange color shows the acoustic measurement of gas bubbles in the water column during the experiment and reported in [68] Environ Fluid Mech to the source. The BPM predicts a trap depth of 740 m, which is very close to the predicted trap depth of 745 m for the CDOG model. For the case with marine diesel in a stronger current (not shown), the BPM predicts a trap depth of 795 m, which lies between predictions from Johansen et al. [67] of 674 m and Chen and Yapa [27] of 824 m. The trajectories of the bubbles above the separation height agree well with the acoustic measurements, with the bubbles rising to similar heights as the observations and following the main trends of the acoustic data. Differences between the model trajectories and the acoustic signals are similar across all models (BPM, DeepBlow, and CDOG) and can be attributed to using crossflow information from a single ADCP in the vicinity of the measurements which may not represent the true velocities along the bubble trajectories. Hence, we conclude that the BPM in TAMOC is validated to an acceptable level with respect to the available measured data in single-phase and multiphase buoyant jets from laboratory to field scales.

Application: API model inter-comparison test cases
To demonstrate the utility of the TAMOC suite to simulate an accidental oil-well blowout and to compare the SPM and BPM solutions, we apply the models to two test cases defined in Socolofsky et al. [31], which reports the results of a model inter-comparison study for oil-well blowout models. The inter-comparison exercise defined several canonical cases of an accidental oil-well blowout in deep water with the purpose to test the performance of predictions for bubble and droplet size distributions at the release and for the fate of oil and gas in the water column over a 14-day simulation. Here, we show how the modules of the TAMOC suite work together to predict oil-well blowout dynamics and compare the TAMOC predictions to those by the other models included in the inter-comparison study.

Test cases
In this paper, we simulate Cases 1 and 2 of the model inter-comparison test matrix, which involve a release of 20,000 bbl/day of oil. The oil is released at a temperature of 150 C, a Environ Fluid Mech water depth of 2000 m, and a gas-to-oil ratio (GOR) of 2000 Nft 3 =bbl. The release is assumed to be from a broken vertical riser, with an outlet diameter of 0.3 m. In both cases, the ambient velocity profile is assumed uniform over the depth at a velocity of 0.05 m/s. Ambient vertical profiles of salinity, temperature, and oxygen were taken from the worldaverage profile in Sarmiento and Gruber [69] as reported in Socolofsky et al. [31].
The complete composition of the oil and gas mixture is given in Socolofsky et al. [31] and includes carbon dioxide, nitrogen, methane, ethane, propane, i-butane, n-butane, ipentane, n-pentane, hexane and one pseudo-component labeled as C7?. We apply the DPM in TAMOC to predict the behavior of the gas and liquid phases of this mixture. The release flow rate of 20,000 bbl/day with a GOR of 2000 is a statement of the equilibrium conditions at atmospheric pressure. These rates are converted to the initial conditions at the source by performing an equilibrium flash calculation in the DPM to estimate the mass fractions of each component in the gas and liquid phases at the release temperature and pressure. TAMOC predicts mass flow rates for gas of 7.53 kg/s and for liquid of 31.11 kg/s at the release. These values agree well with those reported in Socolofsky et al. [31] using the MultiFlash TM software package of 7.4 kg/s of gas and 34.5 kg/s of liquid petroleum fluid at the release. Under these pressure and temperature conditions, 42% of the released mass of methane and light hydrocarbons (C1-C4) is initially present in the petroleum liquid phase; hence, dissolution from both phases should be tracked in the model. The DPM in TAMOC also predicts the densities of petroleum gas and liquid at the release conditions, yielding 136 kg/m 3 for gas and 600 kg/m 3 for liquid. The calculated values reported in Socolofsky et al. [31] were 132 kg/m 3 for gas and 599 kg/m 3 for liquid, which also agree favorably with the TAMOC DPM. As these fluids rapidly cool and dissolve into the ambient receiving water, thermodynamic properties and heat and mass transfer rates are continuously computed by the DPM. Thus, for these highly dynamic and real fluids, the complete equations-of-state in the DPM are essential for accurate tracking of the oil and gas evolution.
The remaining initial conditions for the model simulations are the bubble and droplet size distributions at the source. For the model inter-comparison study, blowout simulations for each model were conducted with the same median gas and oil droplet sizes in order to focus model differences on oil and gas dynamics processes rather than differences in particle size distributions; hence, we impose the median particle size values reported in Socolofsky et al. [31]. To obtain a complete bubble and droplet size distribution from these characteristic values, we use a Rosin-Rammler distribution with spread coefficient equal to 1.8, as recommended by Johansen et al. [70]. We truncate predictions that lie above the maximum stable droplet and bubble sizes, which we compute from equations in Clift et al. [9] for liquid droplets and Grace et al. [71] for gas bubbles. In Fig. 15 we present the distributions used in our TAMOC simulations. Case 1 is for a blowout without any subsurface dispersant injection (SSDI), and the volume median diameter d 50 is 3.8 mm [31]. Case 2 is for a blowout treated by SSDI. Socolofsky et al. [31] assume that SSDI reduces the interfacial tension of the liquid phase petroleum fluids at the release by a factor of 200 compared to Case 1. With this smaller interfacial tension, d 50 ¼ 0:2 mm in Case 2. Socolofsky et al. [31] assume that SSDI does not affect the gas phase, and a median bubble size of 8 mm is given for both cases.
We use these data to run blowout simulations for both the SPM and BPM. The SPM ignores the ambient crossflow, predicts multiple intrusion layers, and assumes that all gas bubbles and oil droplets rise to the surface in a core plume above the wellhead. The BPM includes the effects of currents, predicts only one intrusion layer, and allows gas bubbles and oil droplets to rise out of the plume or enter the intrusion, depending on their dynamic interactions with the plume. The ambient density profile gives N ¼ 0:0013 s À1 in the bottom 500 m of the water column, and the gas and oil together at the release have a kinematic buoyancy flux of B 0 ¼ 0:6 m 4 =s 3 . For a typical gas bubble rise velocity of 0.2 m/s, the empirical equations in Socolofsky et al. [51] give a peel height h P ¼ 650 m for a Type 1* plume, which means that the gas bubbles will be significantly spread out at the first intrusion. Oil droplets with a rise velocity of less than 0.02 m/s would be expected to intrude based on criteria in Chan et al. [10]. The peel height due to currents alone would be h S ¼ 1; 350 m, based on the correlations in Socolofsky and Adams [11], making this a stratification-dominated plume with h P \h S .
An important question to answer is: Which of the two integral models is the most appropriate to simulate the blowout? Because the plume is stratification dominated, the SPM may be acceptable, but because there is a weak current, the real plume will deflect in the downstream direction, breaking the symmetry of the outer plume, making the BPM acceptable as well. The advantage of the SPM is that it can predict multiple intrusion layers, but it cannot predict any liquid oil droplets that might enter the intrusion-which is predicted to occur based on Chan et al. [10]. The advantage of the BPM is that it can predict the intrusion of small oil droplets, but it remains unclear how to initialize secondary plume structures above the first intrusion-all of the entrained fluid enters the intrusion, and the Type 1* plume will have the gas bubbles and oil droplets spread out over a large horizontal region at the top of the plume, making a secondary plume very diffuse. To understand these trade-offs that occur quantitatively for each model, we simulate Cases 1 and 2 with both models and compare their results with each other and with metrics for other models reported in the model inter-comparison study [31].

Discussion
The model inter-comparison study defined several quantitative metrics as a basis for model comparison. Results for these metrics for the SPM and BPM simulations are given in Table 1, along with the ranges spanning the maximum and minimum values predicted for each metric from the other models reported in [31]. We consider these results in the context of the different modeled processes in the discussion below.
In both models, the hot petroleum fluids rapidly cool as the plumes entrain ambient seawater, changing the buoyancy of both the dispersed and continuous phases. Results for both models have similar vertical scales for the heat transfer processes. For the BPM, the entrained seawater rapidly heats to a maximum temperature of nearly 7 C in the first 3 m of plume rise. This temperature increase results in a decrease of the water density from 1037.1 to 1036:4 kg/m 3 , which adds a small amount of buoyancy to the plume. Heat transfer from the dispersed oil and gas continues over the first 25 m along with entrainment of colder ambient fluid. At 50 m above the release, the oil and gas have fully cooled to the temperature of the plume fluid, and the multiphase plume continues to rise until ultimately forming an intrusion layer. This rapid heat transfer at the source is similar to that reported in Anderson et al. [72], who also evaluated the potential for methane hydrate formation in these deepwater plumes. For Case 2, the vertical scales over which heat transfer occurs are similar, but the maximum temperature in the plume reaches 14 C, due to the faster heat transfer from the smaller oil droplets, which have a greater surface-area-to-volume ratio. Hence, heat transfer is a rapid process that has a negligible effect on the final temperature of the seawater and can be neglected as long as the discharged petroleum fluids are initialized with their properties in equilibrium with the ambient seawater temperature.
Both the SPM and BPM predict similar trap depths in Cases 1 and 2 for the first intrusion, which are also within the range predicted in the model inter-comparison study ( Table 1). This is expected since both model formulations have been calibrated to trap depth of single-phase plumes in stratification. The trap depth for the BPM is slightly lower than the SPM in Case 1 because of the extra forced entrainment in the BPM. Both models predict trap depths slightly below the theoretical value from Socolofsky and Adams [51] (1620 m depth), resulting from a combination of the non-linear stratification and the dissolution, which are accounted for in the SPM and BPM and neglected in the theoretical equation. We plot a profile through the BPM results for Case 1 in Fig. 16 to illustrate the intrusion dynamics and subsequent bubble and droplet rise. The weak crossflow bends the plume over in the downstream direction, and the descending intrusion flow overlaps the previous upward-rising plume solution (see the loops in the line depicting the lower edge of the entrained fluid boundary). This overlap occurs because the plume is stratification dominated, so that the currents are too weak to advect the intruding flow away from the upwardrising buoyant jet. For this crossflow, the overlap region occurs over a limited vertical extent (about 50 m out of 300 m total plume rise). Because the BPM predicts a trap depth similar to the SPM, we conclude that this small overlapping flow region, which is neglected in the BPM, has a negligible effect on the dilution or dynamics of the BPM simulation. Hence, even in stratification-dominated conditions, when weak crossflows occur, the BPM gives reliable results for the first intrusion layer.
Comparing the trap depths between Cases 1 and 2 shows that the plumes with droplets treated by SSDI intrude at a slightly lower depth (i.e., they rise higher in the water column). There are two aspects of the models that are affected by the droplet size. First, as explained above, the smaller droplets in Case 2 cool and dissolve faster, which decreases their buoyancy. This effect is included in both models. In both models, this lower buoyancy results in a smaller velocity of the plume fluid, leading to less shear entrainment, which causes it to take longer to reach a neutral buoyancy level, allowing it to rise higher in the water column. Second, the smaller oil droplets stay closer to the plume centerline in Case 2 than the larger droplets of Case 1, which move upstream, away from the plume centerline. This effect is only included in the BPM as the SPM neglects crossflow. In Case 1, as the oil droplet move off the centerline, their buoyant force is removed from the plume solution by j. In Case 2, the smaller oil droplets remain closer to the plume centerline, producing a higher buoyant force in the model equations. This effect will accelerate the plume, leading to higher net entrainment, but also allowing the plume to achieve a greater peel height (shallower depth) in Case 2 than Case 1. By these processes, both the SPM and BPM predict that treated plumes will trap at slightly lower depths than untreated spills.
The effect of the different entrainment algorithms in the models is more evident in the volume flow rate at the first intrusion. Both models predict values within the range of values in the model inter-comparison study (refer to Table 1), though consistently on the Fig. 16 Cross-sectional view of the bent plume model and the particle trajectories after 14 days without dispersant treatment (Case 1) Environ Fluid Mech lower end of the range in Socolofsky et al. [31]. Also, the BPM predicts higher volume flow rates to the intrusion than the SPM. This results for several reasons. The BPM includes forced entrainment, which is neglected in the SPM. The SPM exchanges fluid between the inner plume and downdraught plume, reducing the total volume flux to the intrusion, a process likewise neglected in the BPM. Also, the volume flow rate in the intrusion for the BPM grows rapidly near the neutral buoyancy level, and the value predicted by the model is sensitive to the point where the simulation is terminated. Here, we terminate the BPM at the first neutral buoyancy level in the intrusion. For Case 2, these differences are also compounded by the extra buoyancy provided by the smaller droplets that remain near the plume centerline in the BPM. As a result of these complexities for predicting intrusion layer flow rate, integral plume models are much more accurate for predicting trap depth than dilution at the start of the intrusion layer.
From a response perspective, one of the most important aspects of the simulations is the fate of the petroleum fluids in the water column. Referring back to Fig. 16, the trap depth of the first intrusion is close to the seafloor, and most of the dispersed phase transport will either be in the formation of subsequent intrusion layers (SPM) or Lagrangian particle transport in the water column (BPM passing particles to LPM). The observations during the Deepwater Horizon accident suggest that at most two subsurface intrusion layers resulted from the buoyant plume of oil and gas released at the seafloor [12,73]. Also, during Deepwater Horizon, the oil did not all surface directly above the wellhead, but rather surfaced with an offset in the down-current direction [74,75]. Hence, for tracking oil in the water column, the SPM is not the best tool. As depicted in Fig. 16 for Case 1, the BPM predicts that some oil droplets enter the intrusion and that most of the released oil and gas rises out of the top of the plume along Lagrangian pathways tracked by the LPM. Table 1 reports the mass flow rate of liquid petroleum entering the intrusion layer and reaching the sea surface. In both Cases 1 and 2, the particles initially released as gas bubbles eventually condense to liquid phase, and because these particles initially contain C7?, none of the particles released at the discharge fully dissolve in the water column. In Case 1, the SPM predicts 78% of the released mass surfaces above the wellhead, and in Case 2, 73% reaches the surface. Because the SPM does not allow the dispersed phase to intrude, the missing mass was dissolved in the intrusions, and the smaller droplets of Case 2 give more dissolution. For the BPM in Case 1, most of the released fluid also reaches the sea surface, with the BPM predicting 3.6% entering the intrusion layer at the near field plume termination. After 14 days of simulation, the BPM predicts 75% reaching the sea surface and 25% becoming dissolved in the water column. In Case 2, the smaller droplets lead to more oil entering the intrusion and a greater fraction being dissolved. In the BPM predictions, 69% of the mass initially enters the intrusion, and only 31% continues toward the sea surface near the wellhead. After 14 days, some of the droplets initially in the intrusion have also surfaced so that 37% of the released fluids are reaching the sea surface (though some droplets surface very far downstream), with the remaining mass either dissolved (30%) or suspended (33%) in the water column. Each of the values from these models fall within the ranges predicted by other models in the model inter-comparison study [31]. The important difference between the SPM and BPM predictions are the surfacing locations: the SPM oil surfaces directly over the wellhead, and some of the droplets in the BPM simulations surface many kilometers downstream (Fig. 16).
Because there are no real data for these inter-comparison test cases, the results only serve to compare TAMOC to other response models used in industry. Based on Table 1, both models (SPM and BPM) give trap depths and intrusion layer flow rates within the range of other models. The main weaknesses of the SPM is that it neglects the crossflow and does not allow small oil droplets to intrude. From a response perspective, this is a major disadvantage. However, the BPM is not strictly valid in the typically stratificationdominated conditions of a deepwater blowout where the intrusion low will overlap the upward-rising portion of the plume. Because both models give similar results for the trap depth, we may conclude that the limitations on the BPM are relatively small. Under this assumption, the BPM can be used in all deepwater blowout conditions and is likely to yield acceptable predictions for the pathways of gas bubbles and oil droplets, limited in real spill scenarios primarily by the possibly unknown composition of the spilled petroleum fluids, their mass flow rates, and the unsteady, non-uniform currents in the intermediate field of the blowout. These uncertainties on boundary and initial conditions plague all oil spill response models, and the models such as those contained in the TAMOC suite yield a satisfying level of performance compared to measured validation cases over a wide range of multiphase buoyant jets.

Summary and conclusions
In this paper, we present a comprehensive description and validation of the models in the Texas A&M Oilspill (Outfall) Calculator (TAMOC), which includes modules for the dispersed-phase equations of state in the Discrete Particle Model (DPM), particle tracking in the Lagrangian Particle Model (LPM) and multiphase integral plume models based on an Eulerian double-plume model in the Stratified Plume Model (SPM) and on a Lagrangian framework for a crossflow dominated plume in the Bent Plume Model (BPM). Previous papers discuss the validation of the DPM [29] and the SPM [8], and here we focus on validation of the BPM together with a comparison between the BPM and SPM for predicting deepwater accidental oil-well blowouts.
We created the BPM as a synthesis of different modeling approaches for other integral models in the literature. The shear entrainment algorithm follows the equations in Jirka [4] for transition from an accelerating to decelerating plume, and the Lagrangian element solution follows the algorithm in Lee and Cheung [3], including their maximum entrainment hypothesis that the net entrainment should be the maximum of either the shear or forced entrainment. Our algorithms for dispersed phase particle tracking in the plume are new and have been described in detail in the Methods section. Altogether, the complete model has one calibration parameter j that reduces the effective buoyant force of dispersed phase particles as they approach the edge of the Lagrangian plume element. All other calibration parameters (notably, the entrainment coefficients) are kept at their literature values for single-phase plumes. Our numerical solution is also unique, employing an adaptive step-size solver for stiff ODEs.
To validate the BPM, we considered data for single-and multiphase buoyant jets from laboratory to field scales. The model performance in single-phase buoyant jets is very good, and the adaptive step-size solver preserves high accuracy from the early stages of the jet discharge. In unstratified ambients, this allows the BPM to perform better than other models using the same basic equations; however, in stratified environments, the better accuracy of the entrainment leads to a slightly lower plume rise than observed in experiments when using entrainment coefficients calibrated for use by other models. This performance is still very good, and changes in the entrainment coefficient values to match stratified conditions will negatively impact performance in unstratified cases. In any case, the BPM predicts the intrusion depth in stratification within about 10% relative error.
In multiphase flows, we validated the BPM to laboratory experiments of bubble plumes and bubbly jets and to field-scale simulated oil-well blowouts. Comparing the model predictions for the separation height of gas bubbles on the upstream side of a bubble plume in uniform crossflow, we calibrated a quartic expression for j. Using this relation, we obtain good predictions for the trajectories of gas bubbles and entrained plume fluid for bubble plumes in crossflow. This expression was validated to similar experiments for bubbly and oily jets, for a bubble plume in stratified crossflow, and for the DeepSpill field data. The DeepSpill data further provided a test for the LPM since much of the gas bubble and oil droplet transport for that experiment was as Lagrangian particles advecting in the currents above the deep intrusion layer.
In the field, both stratification and crossflow are ubiquitous, and for multiphase plumes where separation can be expected between the entrained ambient fluid and the dispersed phase particles, both the SPM and BPM have advantages and disadvantages. To illustrate these differences, we simulated two test cases from a model inter-comparison study for accidental oil-well blowouts. In weak crossflow, the simulation results confirmed that the BPM will double-back on itself, and in the interaction region between the upward-flowing initial plume and downward-flowing intrusion flow, the model is not self-consistent-this interaction is neglected. Yet, this shortfall has a negligible impact on the trap depth predicted by the model, and the great advantages of the BPM to consider the entrainment from the crossflow and the possibility for small oil droplets to enter the intrusion make it the preferred model over the SPM for field-scale response simulations.
Indeed, during response to a subsea oil spill, the most important aspect of the model is to predict the fate of the spilled petroleum fluids in the water column. The DeepSpill experiments, field observations during Deepwater Horizon, and the model inter-comparison study test cases all point to the fact that the integral plume stage of a blowout is likely to ascend on the order of 250-450 m above the seafloor. Above this region, the plume is very likely to spread out, and subsequent intrusion layers will be damped by the increasing currents with height above the bottom. This means that for oil-well blowouts in deep water (h [ 1500 m), the Lagrangian particle tracking phase of the transport may be the dominant aspect of the fate of petroleum fluids in the water column. In these cases, it is the combination of the particle rise velocity and the ambient currents that controls where the particles will surface, and this depends most importantly on the oil droplet sizes. Smaller droplets rise slower, transporting them farther downstream, and allowing more time and surface area per unit volume for dissolution and biodegradation, than larger droplets, resulting in a greater fraction of their fate to be in the water column. Considering the validations conducted here, we conclude that the models in the TAMOC suite perform very well in their trajectory modeling and will accurately predict dispersed phase fate when the size distributions and compositions of the released bubbles and droplets are known. and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.